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An  adaptive  finite  element  procedure  is  developed  for  the  transient  analysis  of 
nonlinear  shells.  The  scheme  is  an  h-method  which  employs  fission  and  fusion. 
Criteria  based  on  incremental  work  and  deviation  of  the  bilinear  finite  element 
approximation  to  the  shell  from  a  Kirchhoff -Love  surface  are  used  as  criteria  for 
adaptivity.  The  example  problems  show  that  the  adaptive  schemes  are  capable  of 
achieving  substantial  improvements  in  accuracy  for  a  given  computational  effort. 

They  include  both  material  and  geometric  nonlinearities  and  local  and  global 
buckling. 

In  order  to  formulate  an  r-adaptive  method,  the  conservation  laws,  the 
constitutive  equations,  and  the  equation  of  state  for  path-dependent  materials  are 
formulated  for  an  arbitrary  Lagrangian-Eulerian  description.  Both  geometrical  and 
.  material  nonlinearities  are  included  in  this  setting.^  A  Petrov-Galerkin  method  is^, 
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19.  ABSTRACT  ( continued) 


developed  for  the  stress  update  so  that  the  history  dependence  and  the  resulting 
convective  term  on  the  stress  tensor  can  be  treated.  A  collocation-weighted 
residual  scheme  is  also  developed.  In  addition,  the  tangent  stiffness  matrix  for 
the  equilibrium  equation  is  derived  from  the  principle  of  virtual  work.  Various 
methods  for  solving  the  finite  element  equations  are  presented,  and  several 
numerical  examples  are  analyzed  to  examine  some  features  of  the  proposed  methods. 
The  first  are  some  elastic -plastic  wave  propagation  problems  which  serve  to  check 
the  correctness  of  the  numerical  scheme.  The  second  is  a  flexural  problem,  the 
response  of  which  is  dominated  by  the  formation  of  hinge  lines .  The  adaptive  mesh 
technique  enables  this  problem  to  be  solved  with  a  much  coarser  mesh. 
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1 .  INTRODUCTION 

The  nonlinear  transient  analysis  of  structures  is  a  particularly 
promising  field  for  adaptive  procedures,  because,  among  the  various  classes 
of  structural  finite  element  applications,  it  is  computationally  the  most 
demanding.  Furthermore,  it  is  the  class  of  applications  in  which  a  priori 
selection  of  an  appropriately  refined  mesh  is  most  difficult,  since  the 
areas  of  the  mesh  which  need  to  be  refined  depend  on  the  evolution  of  the 
response,  which  cannot  be  foreseen  by  the  analyst.  Thus,  while  expert 
systems  may  prove  to  be  quite  effective  in  guiding  a  user  to  design 
appropriate  meshes  for  linear-elastic,  static  problems,  it  is  doubtful  that 
this  could  be  done  in  a  typical  nonlinear  transient  problem,  such  as  the 
simulation  of  a  high-energy  disposition  on  a  missile  nose  or  a  front-end 
crash  of  an  automobile.  In  this  type  of  analysis,  the  computational  power 
must  be  focused  on  those  parts  of  the  mesh  which  undergo  the  most  severe 
deformation,  such  as  hinging  and  wrinkling,  and  the  sites  of  such 
deformations  are  not  determinable  a  priori.  Furthermore,  it  is  desirable  to 
start  various  types  of  simulations,  such  as  a  frontal  and  side  crash,  with 
the  same  mesh  and  let  the  response  dictate  any  refinement. 

While  nonlinear  transient  analysis  is  one  of  the  most  promising  areas 
for  adaptive  procedures,  it  is  also  the  most  challenging.  Perusal  of  the 
reviews  of  the  adaptive  field  recently  written  by  Noor  and  Babuska  (1987) 
and  Oden  and  Demkowicz  (1988)  reveal  that  the  bulk  of  the  theoretical  work 
has  been  devoted  to  determining  local  error  estimates  for  linear  static 
problems;  these  estimates  are  used  to  select  the  elements  or  subdomains  to 
be  refined.  These  error  estimates  have  evolved  into  two  main  types: 

1.  residual  error  criteria  based  on  the  magnitude  of  the  residual 
in  the  governing  equations; 
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2.  error  indicators  based  on  interpolation  and  extension  methods. 

A  difficulty  in  applying  these  methods  to  shells  is  that  in  the  most 
effective  elements  for  shell  analysis,  namely  bilinear  quadrilateral 
elements  such  as  described  in  Belytschko  et  al  (1984)  and  in  Hallquist  and 
Benson  (1986),  even  the  shape  of  the  shell  is  not  adequately  represented. 

In  other  words,  while  the  residual  for  the  bilinear  description  of  the  shell 
may  indicate  a  small  residual,  the  errors  may  be  quite  large  due  to 
discrepancies  between  the  configuration  of  the  Kirchhoff -Love  shell  surface 
and  the  bilinear  finite  element  representation. 

However,  this  drawback  also  provides  an  opportunity,  for  in  fact  it  is 
in  the  regions  of  maximal  deviation  between  the  bilinear  representation  and 
the  shell  surface  that  the  finite  element  mesh  is  most  inadequate.  Since  an 
average  normal  to  the  shell  surface  can  be  estimated  at  all  times,  the 
deviation  of  the  bilinear  representation  and  a  more  accurate  approximation 
to  the  shell  surface  can  be  computed  and  used  to  indicate  where  mesh 
refinement  will  prove  useful. 

Another  aspect  of  the  approach  taken  here  is  that  we  have  not 
endeavored  to  obtain  a  certain  level  of  accuracy  by  the  refinement.  This 
choice  was  based  on  two  reasons:  (1)  it  is  impossible  at  this  time,  with 
the  available  mathematical  tools  for  error  estimation,  to  estimate  the  local 
error  in  a  nonlinear  transient  solution;  (2)  in  most  computer  systems,  the 
fast  memory  allocated  to  a  run  must  be  set  at  the  beginning  of  the  run. 

Therefore,  the  philosophy  of  the  adaptive  process  described  here  is  to 
obtain  the  most  accuracy  for  a  given  set  of  computational  resources.  As 
will  be  seen  in  the  examples,  this  philosophy  is  quite  effective.  By  using 
an  adaptive  mesh,  it  is  generally  possible  to  obtain  a  solution  of 
comparable  accuracy  with  half  the  total  number  of  elements,  and  hence,  half 
the  computational  resources  in  an  explicit  program,  as  with  a  fixed  mesh. 
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2.  FINITE  ELEMENT  FORMULATION 


The  shape  of  the  midsurface  is  described  in  this  finite  element  model 


x£  “  Nj  (£)  x^j  (2.1) 


where  x^  are  the  coordinates  of  node  I  and  N^(£)  are  the  bilinear 

isoparametric  shape  functions.  Lowercase  subscripts  designate  Cartesian 
components,  and  uppercase  subscripts  designate  node  numbers;  repeated 
indices  are  summed  over  their  range,  4  for  uppercase,  3  for  lowercase. 

The  shape  functions  N^(£)  are  functions  of  the  reference  variables  £., 

i  -  1,2,  also  written  as  $2”  ^ >  an<*  they  are  given  by 


NI  “  4  (1  +  *1°  (1  +  1 


(2.2) 


where  and  are  the  coordinates  of  the  nodes  at  the  corners  of  the 

reference  domain  defined  by  -1  <  £  <  1,  -1  <  tf  <  1.  Note  that  unless  a  mesh 

of  these  elements  is  quite  refined,  they  provide  a  rather  poor  model  for  the 
curved  surface  of  a  shell.  In  a  regular  mesh  on  a  cylindrical  shell,  these 
elements  are  in  fact  all  flat,  and  any  interaction  between  flexure  and 
membrane  response  only  occurs  at  the  nodes.  Furthermore,  in  a  region  of 
large  curvature,  this  model  can  deteriorate  even  more  severely  with  very 
large  angles  between  adjacent  elements.  However,  before  discussing  this 
further,  the  basic  mechanics  of  this  finite  element  formulation  will  be 
described . 
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In  the  formulation,  two  types  of  coordinates  are  used  in  addition  to 


the  global  Cartesian  coordinates: 

AAA 

1.  for  each  element,  an  element  coordinate  system  (x,y,z)  with 
base  vectors  e^ ,  £2,  and  e^  is  defined  so  that  e^  and  e ^  ate 

tangent  to  the  midsurface  and  rotate  with  the  element; 

2.  for  each  node,  a  triad  b^  is  defined  so  that  it  rotates  with 

the  node,  with  b^  normal  to  the  midsurface  of  the  shell  in  the 
undeformed  configuration. 

Whereas  in  the  original  formulation  of  Belytschko  et  al  (1984)  the 
original  orientation  of  b^  was  arbitrary,  it  is  used  here  to  locate  new 

nodes  created  in  the  adaptive  process  and  therefore  must  initially  be 
approximately  normal  to  the  midsurface  of  the  shell. 

The  deformation  of  the  element  is  governed  by  the  Mindlin-Reissner 
hypothesis,  which  allows  transverse  shear  but  requires  the  normal  to  remain 
straight,  so  the  velocity  of  a  generic  point  in  the  shell  is  given  in  terms 

of  the  velocity  of  the  midsurface  v?  and  the  angular  velocity  by 

A 

vi  -  V®  -  z  (e3  x  (2.3) 


A 

where  z,  by  the  definition  of  the  element  coordinates,  is  the  distance  of  a 
point  from  the  midsurface. 

The  velocity  field  is  given  by 


Since  the  element  uses  one-point  quadrature,  the  strains  are  evaluated 
single  point,  the  origin  of  the  reference  plane  which  is  given  by 
The  velocity  strains  at  this  point  are  given  by 
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and  n  is  the  normal  to  the  shell. 


The  nodal  forces  are  computed  from  the  stresses  by  one-point 


quadrature,  which  yields 


where 
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and  h  is  the  thickness. 

The  one-point  quadrature  element  is  rank  deficient,  so  it  is  associated 
with  spurious  singular  modes,  as  described  in  Belytschko  et  al  (1984). 

Their  control  is  also  described  therein. 

The  incremental  work  is  computed  in  each  element  by 


r  ,,n+h.T  .  n  L  n+i,  ... 
-  J  At  (d  )  (a  +  a  )  dV 


(2.10a) 


where 


d  -  [d  ,  d  ,  d  ,  d  ,  d  ] 
1  xx  ’  yy  xy ’  xz ’  yz  ‘ 


(2.10b) 


a  —  [a  ,  a  ,  a  ,  a  ,  a  1 
-  1  xx'  yy  xy’  xz’  yzJ 


(2.10c) 


and  At  is  the  time  increment;  superscripts  indicate  the  time  step.  This 
quantity  is  used  to  check  stability  and  as  a  criterion  for  mesh  adaptation 
in  some  of  the  studies. 


3.  FUSION  AND  FISSION  ADAPTIVITY 


The  type  of  adaptivity  which  has  been  adopted  here  is  an  h-type,  where 
the  mesh  is  selectively  refined  in  parts  of  the  domain  during  the  evolution 
of  the  solution.  In  addition,  the  refined  elements  are  fused  when  they  are 
no  longer  needed,  so  that  computational  power  is  not  wasted  on  those  parts 
of  the  domain  which  no  longer  undergo  a  changing  deformation  pattern.  The 
motivation  for  including  the  fusion  process  is  that  in  transient  nonlinear 
problems,  certain  parts  of  the  domain  in  effect  "freeze",  so  that  coarse 
meshes  can  capture  their  behavior  effectively. 

The  adaptive  process  consists  of  fission,  in  which  an  element  is  split 
into  four,  and  fusion,  in  which  a  group  of  four  elements  is  combined  into 
one.  These  processes  are  illustrated  in  Fig.  1.  For  purposes  of 
organization,  any  group  of  four  elements  which  is  created  by  fission  is 
called  a  molecule. 

There  are  three  aspects  to  the  implementation  of  fission- fusion 


adaptivity: 


1.  criteria  for  fission  and  fusion;  the  evaluation  of  these 
criteria  is  called  a  judgment; 

2.  the  initial  conditions  for  element  and  nodal  variables  at  the 
nodes  and  elements  which  are  created  by  fission; 

3.  the  initial  conditions  for  element  variables  of  elements 
created  by  fusion. 


Fission-Fusion  Criteria 


Two  criteria  have  been  adopted  for  making  judgments  on  fission- fusion 


1.  an  incremental  internal  work  criterion; 


2.  a  discontinuity  criterion  based  on  the  increase  in  the  angle 
between  two  adjacent  elements. 

According  to  the  incremental  internal  work  criterion,  the  elements 
which  are  fissioned  are  those  which  sustain  the  most  work.  Because  this 
variable  usually  has  an  oscillatory  character  in  an  explicit  solution  of  a 
transient  problem,  the  judgment  is  made  on  the  basis  of  the  total 
incremental  work  done  over  the  last  five  time  steps.  For  the  purpose  of 
comparing  elements  of  different  sizes,  the  total  incremental  energy  in  a 
molecule  is  used  as  the  criterion.  Thus,  fusion  is  indicated  whenever  the 
incremental  work  in  a  group  of  four  elements  which  has  been  created  by  a 
previous  fission  is  smaller  than  the  incremental  energy  in  other  molecules. 

Even  with  the  filtering  that  is  brought  about  by  taking  the  incremental 
work  over  five  steps,  the  incremental  work  criterion  can  lead  to  an 
oscillatory  pattern  of  fission  followed  by  fusion  in  many  molecules  in  a 
transient  process.  Therefore,  a  time  delay  has  been  included  which  prevents 
fission  or  fusion  unless  it  is  indicated  by  two  consecutive  judgments.  This 
type  of  retardation  of  the  adaptive  process  appears  to  be  needed  in  explicit 
treatments  of  nonlinear  structural  dynamics  with  adaptive  meshes  if 
excessive  "churning"  between  fission  and  fusion  is  to  be  avoided. 

The  second  criterion  we  have  studied  is  based  on  the  change  in  angle 
between  two  elements.  The  basis  for  this  criterion  is  that  one  of  the 
largest  sources  of  errors  in  this  finite  element  procedure  is  the  inability 
of  the  piecewise  bilinear  elements  to  capture  the  correct  shape  and  moment- 
curvature  interaction  of  the  shell  as  the  deformation  localizes.  Severe 
deformation  in  shells  is  usually  associated  with  large  curvatures:  since  the 
bilinear  element  cannot  represent  large  curvatures  directly,  it  is 
associated  with  severe  "kinking"  between  elements,  which  can  be  detected  by 
monitoring  the  angle  between  elements.  For  those  elements  which  satisfy  the 


angle  criterion  for  fission,  both  elements  on  each  side  of  the  line  are 
subdivided  into  four  elements . 


An  advantage  of  the  angle  criterion  over  the  incremental  work  criterion 
is  that  it  can  be  applied  to  more  than  one  level  of  fission- fusion  without 
need  of  additional  parameters.  The  incremental  work  criterion,  if  it  is  to 
be  used  for  two  levels  of  fission-fusion,  requires  the  specification  of  a 
ratio  at  which  the  second  level  of  fission  is  initiated. 

Since  it  is  difficult  to  relate  any  of  these  criteria  to  the  ultimate 
accuracy  of  a  solution,  one  technique  we  have  frequently  used  is  to  simply 
specify  the  maximum  number  of  elements  and  use  the  criteria  to  select  where 
those  elements  are  placed.  In  this  procedure,  we  start  with  a  uniform  mesh 
which  contains  a  fraction  of  the  maximum  number  of  elements  allowed.  After 
five  time  steps,  the  elements  are  fissioned  in  decreasing  order  of  the 
amplitudes  of  their  error  indicators  until  the  maximum  number  of  elements  is 
obtained. 

Fission  Process 

When  a  molecule  is  fissioned  next  to  an  unfissioned  molecule,  as  shown 
in  Fig.  2,  nodes  are  created  adjacent  to  unsplit  sides,  so  they  cannot  be 
handled  by  the  usual  equations  of  motion.  In  order  to  correctly  handle 
compatibility,  these  nodes  must  be  treated  as  "slave”  nodes  which  are  driven 
by  the  adjacent  "master"  nodes.  In  addition,  in  order  to  introduce  a  good 
representation  of  the  shell  as  quickly  as  possible,  it  is  useful  to  use  the 
nodal  vectors  b^  for  an  approximation  of  the  curved  Kirchhoff -Love  surface 

on  which  the  new  nodes  are  placed. 

The  procedure  for  setting  the  initial  configuration  of  the  nodes  is  as 
follows.  The  surface  is  approximated  oy 


where  Hj(£)  are  the  Hermite  interpolants ,  so  that  Hj  ^(Cj)  -  5jj  and  4*^  is 

the  slope  of  the  Kirchhoff  surface  relative  to  the  bilinear  approximation, 
which  is  obtained  by 

^IJ  "  '^-3  ’  — I HlJ I  (3.3) 

The  initial  velocities  of  the  nodes  are  obtained  from  the  bilinear 


interpolation  (2.1).  The  initial  element  variables  for  the  elements  are 
taken  from  the  parent  element.  The  mass  matrix  is  reassembled  after 
fission.  Nodes  which  are  formed  at  sides  which  are  continuous  sides  of  the 
adjacent  elements  are  considered  slave  nodes.  All  other  new  nodes  are 
master  nodes.  Thus,  an  interior  node  is  always  a  master  node,  but  if  only  a 
single,  isolated  molecule  undergoes  fission,  all  other  new  nodes  are  slave 


In  the  fusion  process,  no  new  nodes  are  created;  the  velocities  and 
displacements  at  the  nodes  which  remain  are  assumed  to  be  continuous  during 
fusion.  The  number  of  elements  is  reduced  from  four  to  one;  the  historical 
state  variables  (stress  components  and  yield)  are  taken  to  be  the  area- 
weighted  average  of  the  parent  element  stresses.  State  variables  such  as 
the  yield  stress  are  adjusted  so  they  remain  consistent.  For  example,  if 
the  majority  of  the  elements  is  plastic,  the  yield  stress  is  adjusted  so 
that  the  fused  element  is  also  yielding. 


l-’it 


4 .  NUMERICAL  examples 


All  numerical  examples  described  in  this  section  were  performed  on  a 
Harris  800  computer  in  single  precision.  A  single-precision  word  consists 
of  11  significant  digits  (base  10)  on  this  computer. 

Since  closed- form  solutions  are  not  available  for  nonlinear  transient 
problems,  two  types  of  comparisons  are  used  for  the  adaptive  solutions: 

1.  numerical  results  obtained  by  finer  meshes; 

2.  experimental  results. 

The  first  example  concerns  a  clamped  beam  which  is  impulsively  loaded 
over  the  center  portion  as  shown  in  Fig.  3.  Using  symmetry,  half  of  the 
beam  is  modelled  by  m  x  n  quadrilateral  plate  elements,  with  m  elements 
across  the  1.2  in.  width  and  n  elements  over  the  5  in.  half-span.  The  x  and 
z  components  of  the  translations  and  rotations  about  the  x  and  z  axes  were 
constrained. 

Figure  4  shows  the  midspan  deflection  obtained  by  two  fixed  meshes  and 
an  8-  to  10-element  adaptive  mesh.  As  can  be  seen,  the  adaptive  mesh  is 
quite  close  to  the  20-element  fixed  mesh  for  the  first  0.5  msec,  and  it 
matches  the  maximum  displacement  quite  well.  Subsequently,  it  diverges 
somewhat  from  the  fine -mesh  solution. 

Figure  5  shows  the  pattern  of  mesh  adaptivity.  The  first  elements  to 
be  fissioned  are  those  beneath  the  impulsive  load;  the  location  of  the 
fissioned  elements  then  moves  back  and  forth  between  the  center  and  the 
support,  like  the  hinge  in  the  rigid-plastic  solution,  and  finallv  fixes 
itself  at  the  clamped  wall. 

The  profiles  of  the  beam  obtained  by  a  fine  fixed  mesh  and  the  adaptive 
mesh  are  compared  in  Fig.  6.  As  can  be  seen,  the  profiles  compare  t  ..te 
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well,  except  at  the  final  time,  0.49  msec,  when  the  node  at  the  adaptive 
mesh  at  x  -  3.0  in.  deviates  markedly. 

This  beam  was  resolved  with  a  2  x  20  fixed  mesh  and  a  16 -element 
adaptive  mesh.  The  midspan  deflection  is  shown  in  Fig.  7.  In  this  case, 
the  adaptive  mesh  corresponds  very  closely  with  the  fixed  mesh,  even  though 
it  required  only  40%  of  the  elements .  The  potential  savings  in 
computational  resources  is  even  greater  because,  in  the  adaptive  mesh,  half 
of  the  elements  could  employ  a  time  step  twice  as  large  as  that  used  in  the 
fixed  mesh. 

A  more  complex  example  for  the  adaptive  mesh  is  provided  by  the 
cylindrical  panel  problem  shown  in  Fig.  8.  An  initial  velocity  of  5650 
in/sec  is  applied  to  the  3.08  in.  x  10.205  area  indicated  in  Fig.  8.  The 
panel  is  simply  supported  at  its  ends  and  clamped  at  the  sides.  An  Ilyushin 
plasticity  model,  which  is  expressed  in  terms  of  the  resultant  moments  and 
membrane  forces,  and  ^ ,  is  used  in  the  computation. 

Two  adaptive  meshes  were  used  in  the  computation:  a  96-element  mesh 
based  on  a  4  x  8  mesh  of  molecules;  a  218 -element  mesh  based  on  an  8  x  16 
mesh  of  molecules.  The  results  are  compared  to  uniform  fixed  meshes  with 
32,  96,  128,  218,  and  512  elements. 

The  displacement  time  histories  for  the  coarse  adaptive  mesh  are 
compared  to  three  of  the  fixed-mesh  results  in  Figs.  9  and  10  at  points  A 
(z  -  -6.28  in.)  and  B  (z  -  -9.42  in.)  which  are  indicated  in  Fig.  8. 
Remarkably,  the  53-eleraent  adaptive  result  almost  coincides  with  the  128- 
element  fixed-mesh  result  for  the  first  0.4  msec.  Subsequently,  the  results 
zi  the  two  meshes  deviate  somewhat,  and  the  adaptive  results  become  slightlv 
r:ugh,  which  is  caused  by  excessive  churning  of  the  f i ss ion- fus ion  process 


Note  that  the  fixed-mesh  results  with  fewer  elements  deviate  substantially 
from  the  finest  fixed-mesh  results. 

The  displacements  for  the  fine  adaptive  mesh  are  compared  to  the  fixed 
meshes  in  Figs.  11  and  12  at  points  A  and  B,  respectively.  Here  the  218- 
element  adaptive  mesh  corresponds  quite  closely  with  the  512  uniform  fixed 
mesh  and  exhibits  marked  improvement  over  a  200-element  fixed  mesh. 

Deformed  mesh  plots  for  the  finer  adaptive  mesh  are  shown  at  various 
times  in  Fig.  13.  Here  the  incremental  energy  was  used  for  the  fission- 
fusion  criterion.  It  can  be  seen  that  after  0.0125  msec,  the  crown  settles 
downward  like  a  plateau  and  the  fissioning  process  migrates  laterally 
towards  the  line  where  the  curvature  is  maximum.  During  this  time,  the 
crown  moves  down  in  a  frozen  plateau-like  state.  After  that,  the  crown 
develops  a  convex  curvature  when  viewed  from  above,  and  the  elements  in  the 
crown  are  again  fissioned.  The  end  of  the  simulation  again  exhibits 
churning  of  fission- fusion ,  which  is  a  tendency  that  needs  to  be  fixed:  it 
is  probably  due  to  the  fact  that  incremental  work  is  quite  small  in  the 
later  stages  because  most  of  the  deformation  has  taken  place,  so  the 
incremental  work  in  molecules  is  quite  uniformly  distributed,  allowing  the 
fission- fusion  process  to  be  triggered  by  small  oscillations  in  the 
solution. 

Experimental  results  have  been  obtained  for  this  shell  by  Morino  et  al 
(1971),  who  reported  a  maximum  deflection  of  1.24  in.  at  point  A.  The 
finest  fixed  and  adaptive  meshes  yielded  maximum  deflections  of  1.20  and 
1.17  in.,  respectively. 

The  deformed  profiles  are  shown  in  Figs.  14  and  15  for  the  512  -  e  1  emer.t 
fixed  mesh  and  the  218-element  adaptive  mesh.  The  development  of  a  hinge  - 
like  pattern  at  about  x  -  2.0  in.  and  the  attendant  fission  process  are 
quite  clearly  seen  in  Fig  14 .  Figure  15  shows  a  cross  section  in  *.-.e  v  - 1 
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plane  of  symmetry.  The  fission  process  which  takes  place  while  the  crown  is 
moving  like  a  flat  plateau,  followed  by  the  fission  which  develops  when  the 
crown  curves,  is  quite  clearly  seen. 

The  third  example  is  a  hollow,  cylindrical  column  which  is  subjected  to 
a  compressive  axial  load.  This  problem  is  of  interest  because  it  exhibits 
both  global  and  local  buckling,  the  latter  resulting  in  buckling  of  the 
cross  section.  Numerical  results  and  experimental  results  have  been 
reported  for  this  problem  by  Kennedy  et  al  (1986).  The  problem  parameters 
are  given  in  Table  1. 

The  cylinder  is  loaded  by  prescribing  an  upward  velocity  of  500  in/sec 
to  the  bottom  nodes  of  the  model,  with  the  top  fixed.  To  trigger  the 
lateral  buckling  mode,  an  imperfection  given  by 

Ax  -  0.01  sin 

where  i  is  the  length  of  the  column  and  z  is  the  coordinate  along  the  axis 
of  the  column,  is  added  to  the  x-coordinate  of  all  nodes.  The  pattern  of 
adaptivity  is  shown  in  Fig.  16.  Initially,  the  fission  process  moves  up  and 
down  the  column  similar  to  a  reflected  wave.  The  fission  process  then 
coalesces  at  the  nodes  of  the  lateral  buckling  mode,  where  they  remain, 
except  for  the  fission  which  takes  place  at  the  compressive  buckles  at  the 
top  and  bottom  of  the  column.  The  displacements  of  two  points  in  the 
adaptive  mesh  are  compared  to  corresponding  points  in  a  finer  fixed  mesh  in 
Fig  17.  The  results  show  good  agreement. 


5.  CONCLUSIONS 


An  adaptive  procedure  based  on  an  h-scheme  has  been  developed  for  the 
nonlinear  transient  analysis  of  shells.  This  includes  the  development  of 
suitable  criteria  for  fission  and  fusion  and  the  formulation  of  the  fission 
and  fusion  processes.  Two  criteria  have  been  found  useful  for  the  bilinear 
quadrilateral  elements  commonly  used  in  transient  analysis  by  explicit  time 
integration: 

1.  the  incremental  internal  work  criterion; 

2.  the  relative  angle  criterion,  which  is  a  measure  of  the 
deviation  of  the  bilinear  surface  from  the  Kirchhoff -Love 
surface  associated  with  the  nodal  orientations. 

Furthermore,  to  avoid  excessive  churning  of  the  fission-fusion  process,  time 
delays  had  to  be  incorporated  in  the  judgment  process.  Nevertheless, 
churning  becomes  a  problem  with  the  incremental  work  criterion  in  the  later 
stages  of  impulsively  loaded  problems  when  the  work  on  the  system  decreases. 

The  results  we  have  obtained  show  that  the  adaptive  schemes  are  capable 
of  achieving  substantial  improvements  in  accuracy  for  a  given  computational 
effort.  Generally,  an  adaptive  mesh  is  capable  of  achieving  the  same  level 
of  accuracy  as  a  fixed  mesh  with  less  than  half  of  the  computational 
resources.  The  fission  process  tends  to  take  place  in  the  subdomains  where 
the  maximum  deformation  occurs. 

The  h-adaptive  procedure  is  limited  in  its  ability  to  focus  on  the 
subdomains  of  maximum  deformation  by  the  fact  that  the  molecules  which  are 
subdivided  are  fixed  in  the  reference  configuration.  Therefore,  hinge  lines 
which  occur  at  small  angles  relative  to  the  mesh  lines  are  not  captured 
effectively.  For  this  reason,  an  h-r  adaptive  procedure  which  permits 
motion  of  the  nodes  is  now  under  development.  Ar.  essential  ingredient  or 
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Fusion  and  fission  of  a  molecule  with  the  node  numbering  convention. 


Interface  between  a  fused  and  a  fissioned  molecule  with  curved 
geometry. 
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5000  in/sec;  thickness  -  0.125  in. 


4.  Center-point  deflection  of  the  clamped  beam. 


5.  Undeformed  and  deformed  plots  for  the  clamped  beam  with  a  10-element 
adaptive  mesh. 


Deformed  cross-sectional  profile  of  the  clamped  beam  with  2  x  10  mesh. 
Center-point  deflection  of  the  clamped  beam  with  finer  mesh. 
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13.  Undeformed  and  deformed  plots  for  the  cylindrical  panel  with  16  x  32 
adaptive  mesh. 
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15. 


Deformed  cross-sectional  profiles  of  the  cylindrical  panel  with  16  x  32 
mesh  and  in  an  x-y  plane  passing  through  node  A. 


Deformed  cross-sectional  profiles  of  the  cylindrical  panel  with  16  x  32 
mesh  and  in  a  y-z  plane  passing  through  nodes  A  and  B. 


Undeformed  and  deformed  adaptive  meshes  for  the  cylindrical  column. 
Displacement  time  histories  for  two  points  of  the  cylindrical  column. 
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APPENDIX 


CHAPTER  1 


INTRODUCTION 


The  analysis  of  shells  into  the  failure  domain  or  post-buckling  domain  is 
becoming  of  increasing  importance  in  the  design  of  certain  classes  of  struc¬ 
tures,  particulary  when  the  survivability  or  vulnerability  of  the  structure  is 
a  key  question.  Areas  of  design  in  which  failure  analysis  is  important  are 
the  design  of  defense  structures  and  in  the  analysis  of  re-entry  vehicles.  In 
both  cases  the  loads  to  be  sustained  are  extremely  high  or  the  cost  of  super¬ 
fluous  strength  are  severe,  so  it  is  necessary  to  be  able  to  predict  the 
behavior  of  the  shell  substantially  beyond  the  classical  buckling  load, 
because  buckling  by  itself  does  not  constitute  failure  of  the  structure.  Many 
structures  are  in  fact  designed  to  buckle  and  the  key  question  in  the 
survivability  or  viability  of  such  structures  is  whether  the  final 
displacements  and  deformations  are  within  design  limits.  The  prediction  of 
the  displacements  of  a  shell  in  such  severe  environments  requires  an  analysis 
of  the  shell  structure  into  the  post-buckling  regime. 

The  analysis  of  shell  structures  in  the  post-buckling  regime  still  poses 
formidable  difficulties.  The  major  source  of  these  difficulties  is  the  fact 
that  in  the  post-buckling  regime,  high  strains  and  large  deformations  are 
often  localized  in  small  regions  of  the  shell.  Unless  very  refined  meshes  are 
used  in  the  vicinity  of  these  buckles,  large  errors  occur  in  the  predicted 
deformation  of  the  shell,  which  detracts  severely  from  the  usefulness  of  the 
analysis. 

This  work  seeks  to  remedy  these  difficulties  by  developing  adaptive  mesh 
techniques  for  analyzing  the  post-buckling  behavior  of  shells.  The  proposed 
work  consists  of  two  major  tasks: 


1.  the  development  of  an  effective  technique  for  moving  and  refining  the 
mesh  so  that  the  post-buckling  response  of  shells  can  be  effectively 
analyzed; 

2.  the  development  of  methods  for  integrating  the  stress  tensor  in  time  in 
history-dependent  materials  in  meshes  which  move  relative  to  the 
material . 

This  report  addresses  the  second  question.  In  particular,  some  novel 
methods  have  been  developed  which  can  effectively  treat  moving  meshes  for 
path-dependent  materials.  The  resolution  of  this  difficulty  represents  a 
major  step  toward  the  ultimate  goal. 

These  adaptive  mesh  techniques  in  transient  problems  are  inherently  by 
nature  neither  Lagrangian  nor  Eulerian  but  arbitrary-Eulerian-Lagrangian  (ALE) 
and  we  will  use  this  terminology  here.  This  concept  was  first  proposed  in  Noh 
[1964]  under  the  name  "Coupled  Eulerian-Lagrangian."  Similar  applications  to 
compressible  flow  problems  are  reported  by  Trulio  [1966].  Later  in  HIrt  et 
al.  [1974],  the  ALE  method  Is  employed  in  conjunction  with  an  implicit 
formulation  for  the  solution  of  two-dimensional  flows.  In  that  paper,  the 
calculations  for  each  time  step  are  separated  into  three  phases.  The  first 
phase  consists  of  an  explicit  Lagrangian  calculation  for  the  velocities  and 
specific  internal  energies.  Secondly,  a  Lagrangian  implicit  iteration 
procedure  adjusts  the  state  variables  based  on  the  predicted  variables 
obtained  in  the  first  phase.  This  implicit  procedure  eliminates  the  usual 
Courant-type  numerical  stability  condition.  Finally,  in  the  rezoning  phase, 
the  convective  fluxes  for  the  conservation  equations  are  computed  to  account 
for  the  mesh  motions.  An  extended  version  of  this  computational  technique  to 
three-dimensional  flows  contained  within  arbitrarily  shaped  or  moving 
boundaries  is  reported  in  Pracht  [1975]  and  Stein  et  al.  [1977]. 


The  ALE  method  was  Introduced  into  the  finite  element  method  in  Donea  et 
al.  [1977]  and  Belytschko  and  Kennedy  [1978]  in  response  to  the  need  for 
nonlinear  simulation  techniques  in  nuclear  safety  analyses.  The  advantages  of 
the  ALE  method  in  fluid-structure  interaction  are  apparent  since  the  fluid 
domain  can  be  teated  by  the  ALE  formulation  and  the  structural  domain  can  be 
handled  by  the  usual  Lagrangian  description.  In  these  articles,  the  effort  is 
primarily  directed  toward  inviscid  compressible  fluids,  while  in  Hughes  et  al. 
[1978],  a  finite  element  procedure  for  viscous  incompressible  flows  and  free 
surface  flows  is  presented  in  conjunction  with  a  general  kinematical  theory 
for  the  ALE  description.  A  similar  formulation  has  later  been  reported  in  Liu 
and  Ma  [1982].  The  capability  of  ALE  method  to  handle  an  expanding  gas  bubble 
immersed  in  a  fluid  has  been  demonstrated  in  Belytschko  et  al.  [1982].  The 
mesh  motions  for  this  problem  are  prescribed  according  to  a  simple  control 
algorithm  in  which  the  boundary  nodes  are  considered  Lagrangian  and  the  mesh 
velocities  for  the  intermediate  nodes  are  lineary  interpolated  between  the 
boundary  velocities. 

Recently,  the  application  of  the  ALE  concept  to  contact  problems  between 
flexible  structures  is  proposed  in  Haber  and  Abel  [1983]  in  which  the 
displacement  vector  is  separated  into  the  Lagrangian  and  Eulerian  parts.  The 
slip  compatibility  conditions  are  met  by  making  the  Lagrangian  displacements 
common  to  elements  on  either  side  of  the  interface.  Separate  Eulerian 
displacements  are  associated  with  elements  on  one  side  of  the  interface  to 
model  the  slip  conditions.  This  concept  is  extended  in  Haber  [1984]  for 
quasi-static  solid  mechanics  and  Haber  and  Koh  [1985]. 

The  articles  cited  above  are  mainly  directed  toward  linear  path- 
independent  materials  like  Hookean  solids  and  Newtonian  fluids.  The  stress 
states  for  these  materials  are  solely  determined  by  the  displacement  or 


velocity  fields.  When  the  ALE  method  is  applied  to  materials  with  memory,  the 
state  variables  for  an  element  are  affected  by  the  migration  of  material 
points  which  may  carry  different  stress  and  strain  histories.  This  difficulty 
arises  because  the  adaptive  or  ALE  mesh  does  not  model  the  same  set  of 
material  points  throughout  the  simulation.  A  similar  situation  is  encountered 
in  Derbalian  et  al.  [1978]  when  the  Eulerian  description  is  employed  for 
plastic  forming  analysis.  In  this  paper,  the  difficulty  is  3ide-stepped  by 
interpolating  the  stress  histories  at  each  incremental  step. 

Obviously,  a  consistent  treatment  for  the  transport  of  the  histories  of 
all  path-dependent  quantities  through  the  mesh  is  necessary  for  the  analysis 
of  these  materials.  It  is  therefore  the  purpose  of  the  current  effort  to 
develop  a  general  formulation  and  an  explicit  computational  procedure  for 
nonlinear  adoptive  finite  element  analyses.  Furthermore,  to  provide  a  better 
understanding  of  the  adaptive  methods  in  nonlinear  mechanics,  the 
linearization  procedure  for  the  equilibrium  equation  is  addressed  in  this  work, 
to  examine  some  features  of  the  method. 

The  scope  of  the  present  investigation  is  arranged  as  follows.  In 
Chapter  2,  the  notations  for  material,  spatial,  and  referential  coordinates 
are  introduced.  The  relationship  between  the  material  and  referential  time 
derivatives  is  reviewed  with  the  definition  of  the  relative  velocity  between 
the  material  and  mesh  velocities.  The  balance  laws,  which  include  the 
continuity  equation,  momentum  equation  and  energy  equation,  constitutive 
equations,  and  equation  of  state  In  the  adaptive  description  are  reviewed. 

The  stress-velocity  product  is  defined  to  circumvent  the  difficulty  associated 
with  treating  the  gradient  of  the  stresses  in  the  constitutive  equation.  The 


weak  formulation  and  matrix  equations  are  derived.  Based  on  the  principle  af 
virtual  work,  the  linearization  procedure  in  the  adaptive  description  is 


performed  for  the  Cauchy  stress  and  velocity  gradient.  The  resulting  rate 
form  for  the  virtual  internal  energy  has  a  strong  resemblance  to  the 
Lagrangian  description.  The  tangent  stiffness  matrices  consist  of  two 
parts:  the  first  one  is  identical  to  the  Lagrangian  description  and  the 
second  part  arises  when  the  ALE  description  is  employed.  The  explicit 
expressions  for  these  two  matrices  are  given,  and  the  latter  for  the 
contributions  from  the  ALE  description. 

Numerical  methods  for  the  ALE  equations  are  detailed  in  Chapter  3. 
Techniques  for  nonlinear  convective  effects,  which  characterize  the  ALE 
description,  are  reviewed.  These  include  the  upwind  method,  Petrov-Galerkin, 
and  Taylor-Galerkin  formulations.  The  difficulty  of  treating  the  gradient  of 
stresses  in  the  ALE  description  is  discussed.  The  remedy  of  this  issue  is 
separated  into  two  parts.  The  first  part  includes  the  construction  of  the 
stress-velocity  product  in  conjunction  with  an  artificial  viscosity  technique 
to  achieve  the  streamline  upwind  effect.  The  determination  of  the  artificial 
viscosity  parameter  is  accomplished  by  recourse  to  the  analytic  solutions  for 
a  one-dimensional  constitutive  eqution.  A  collocation  weighted  residual 
formulation  for  the  constitutive  equation  is  presented  in  the  second  part. 
This  procedure  is  formulated  such  that  it  can  handle  any  number  of  quadrature 
points  for  the  family  of  displacement  elements.  An  efficient  two-quadrature 
point  element  suitable  for  nonlinear  calculations  is  discussed.  The  usage  of 
this  element  to  account  for  both  the  geometric  and  material  nonlinearities  is 
described.  An  explicit  time  integration  algorithm  based  on  the  predictor- 
corrector  method  is  presented.  Several  numerical  examples  are  analyzed  to 
demonstrate  the  effectiveness  of  the  present  development. 
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CHAPTER  2 


ARBITRARY  LAGRANG IAN-EULERI AN  FORMULATION 


OF  FIELD  EQUATIONS 


2.1  Review  of  Governing  Eauations  in  Arbitrarv  Lagrangian-Eulerian 


Description 

The  material,  spatial,  and  referential  coordinates  are  denoted  by 
X,  x  and  respectively.  Throughout  this  dissertation,  standard 
indicial  notation  is  adopted;  lower  case  subscripts  denote  the 
components  of  a  tensor  and  repeated  subscripts  imply  a  summation  over 
the  number  of  space  dimensions  (NSD).  A  comma  followed  by  a  subscript 
designates  the  partial  derivative  with  respect  to  the  corresponding 
spatial  variable. 

The  superposed  doc  and  scar  denote  the  time  derivative  with  the 
material  and  referential  coordinates  fixed,  respectively.  The 
relationship  between  these  two  derivatives  are  expressed  here  as 
[Hughes  et  al . ,  1978] 


(')  -  (*)  +  ct(  ),i 


(2.1) 


where  c^  is  the  relative  velocity  between  the  material  (v.-)  and  mesh 
velocitv  (v,), 
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Che  relative  motion  between  the  material  and  the  mesh.  In  particular, 
when  c  is  chosen  to  be  v,  i.e.,  ^  »  x,  the  familiar  Eulerian  time 
derivative  is  obtained. 

The  equations  which  govern  the  continuum  in  the  ALE  description 
are  the  three  conservation  laws  [Donea,  1977] 
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mass :  p  +  cjP^  *  ”Pvi,i 


(2.3) 
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momentum:  pv^  +  pCjV^  j  *  Tjj  (j  +  b^ 


(2.4) 
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energy:  pe  +  pCi«tl  *  ^ijv(i,j)  +  5a  “  9i,i 


(2.5) 


The  Cauchy  stress  may  be  decomposed  into  the  deviacoric  stress  tensor 
s^j  and  che  hydrostatic  pressure  p 

tij  -  s±i  -  ?61;j  (2.6 


which  are  given  by  Che  race  constitutive  equation 
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and  che  equation  of  state 
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respectively.  In  the  above,  p  is  the  density;  is  the  body  force  per 
an it  volume;  e  is  the  internal  energy;  is  the  Kronecker  delta; 
p(p,e)  is  the  function  for  the  equation  of  state;  v(i,j)  is  che  velocity 
strain  tensor 


v(i,j)  *  2  (vi.i  *  vj,lJ 

and  v[i(j]  i3  c^e  spin  tensor 

v[i,j]  "  2  (vi,j  '  vj,i} 


(2.9) 


(2.10) 


a  is  the  internal  heat  generation;  q^  is  the  heat  flux;  and  is 

the  material  response  tensor  which  relates  any  frame-invariant  race  of 
the  Cauchy  stress  [Prager,  1961]  to  Che  velocity  strain.  Both 
geometric  and  material  nonlinearities  are  included  in  the  setting  of 
cq  s •  (  2  • 3—8 ) . 

REMARK  2.1.1:  The  right  hand  sides  of  Eqs.  ( 2. 3-5 , 2. 7-8 )  remain  the 
same  for  all  descriptions  [Liu,  1984]. 

REMARK  2.1.2:  Eqs.  (2. 3-5, 2. 7-8)  are  referred  as  Che  "quasi-Eul erian” 
description  in  Belytschko  et  al.  [1980]  because  these  equations  have  a 
strong  resemblance  to  the  Eulerian  equations.  In  particular,  the 
Eulerian  equations  can  be  readily  obtained  by  choosing  c  »  v,  i.e., 

i  -  *• 


REMARK  2.1.3:  Eq.  (2.7)  is  equivalent  to  Che  following  equations: 
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sij  +  yij k,k  "  ck,ksij 


Cijklv(k,Z)  +  skjv[i,k]  +  skiv[j,k] 


and 


71  jk  "  sijck 


(2.12) 


where  is  che  stress-velocity  product.  In  the  following  finite 

element  computation,  these  two  equations  will  replace  Eq.  (2.7)  in  the 
weak  form;  see  Section  3.3  for  a  discussion. 


The  variational  equations  corresponding  to  the  conservation 

equations,  Eqs.  (2. 3-5, 2. 7-3) ,  are  obtained  by  multiplying  by  the  test 

functions,  dp ,  dvlt  de,  d3ij,  dy^j^  and  dp,  over  the  spatial  domain 

and  employing  the  divergence  theorem  to  imbibe  the  traction  force 

vector  h  on  the  boundary  rx  the  amount  of  heat  3  transmitted 
9 

through  rx* 
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momentum 
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energy 
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constitutive 
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equation  of  state 
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Eq.  (2.13)  and  Eq.  (2.15)  represent  the  control  volume  forms  of 
material  and  energy  conservation,  respectively.  Eq.  (2.14)  is  a 
generalization  of  the  principle  of  virtual  work  to  the  control  volume 
form  with  the  first  integral  brought  in  as  d'Alembert  forces. 

In  finite  element  methods,  the  domain  of  interest  Rx  is  subdivided 
into  elements.  Different  sets  of  shape  functions,  £,  ji3  ,  ^e,  ^s, 

2'7  and  ,  and  corresponding  sets  of  test  functions,  J,  , 

and  are  introduced  to  interpolate  the  velocity,  density, 
internal  energy,  deviatoric  stress,  stress-velocity  product  y  and 
hydrostatic  pressure,  respectively.  Note  that  the  test  functions  and 
shape  functions  for  deviatoric  stresses  are  used  only  in  the 


constitutive  equations.  The  matrix  equations  corresponding  to  Eqs. 
(2.13-18)  are: 
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energy 
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(2.22) 


M^y  =■  N^s 


equation  of  state  HP£  +  Np£  *  u 


(2.23) 


(2.24) 


The  superscript  "T"  denotes  matrix  transpose;  M,  Mp ,  M®,  Ms ,  and  Mp 
are  the  generalized  mass  matrices  for  the  corresponding  variables  in 
Eqs.  (2.19-24),  respectively;  tf,  N*3 ,  Ne,  xj  and  are  the  generalized 
convective  matrices;  Kp  is  the  stiffness  matrix  for  density;  fint  is 
is  the  internal  force  vector;  fexc  is  the  external  load  vector;  £  is 
Che  conductance  vector;  r  is  the  generalized  energy  source  vector;  G  is 
the  divergence  operator  matrix;  D  is  the  generalized  diffusion  matrix 
for  the  deviacoric  stress;  z  and  u  are  the  generalized  deviacoric 
stress  and  pressure  vectors,  respectively.  The  definitions  for  these 
matrices  and  vectors  are  given  in  Appendix  A. 
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REMARK  2.2. 1 :  The  nonlinear  convective  terms,  which  characterize  the 
ALE  method,  inevitably  pose  difficulties.  Recently,  finite  element 
methods  for  non-self-adjoint  systems  (see  Section  3.2  for  a  discussion) 
have  been  developed  which  do  not  suffer  from  crosswind  diffusion  when 
applied  to  the  multi-dimensional  convection-diffusion  problems.  These 
methods  may  be  applied  to  handle  the  convective  terms  in  Eqs. 
(2.19-21,2.23-24). 

REMARK  2.2.2:  All  the  matrices  and  vectors  defined  in  Appendix  A  are 
integrated  over  the  spatial  domain  Rx  which  changes  continuously 
throughout  the  computation. 

REMARK  2.2.3:  The  stress-velocity  product  v  is  scored  at  each  node  as 
a  vector  with  a  dimension  of  (number  of  space  dimensions )x( number  of 
stress  components).  The  diagonal  form  for  is  obtained  by  locating 
the  numerical  integration  points  at  the  nodes. 

REMARK  2.2.4:  A  procedure  for  the  stress  update  equations  (2.22-23)  is 
presented  in  the  next  chapter  to  clarify  the  temporal  integration  for 
path-dependent  materials.  All  the  path-dependent  quantities  are 
updated  analogous  to  Eqs.  (2.22-23).  To  the  author's  knowledge,  it  is 
a  new  approach  to  calculate  the  stress  states  for  path-dependent 
materials  in  the  ALE  description. 
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2.3  Finite  Elemenc  Formulation  by  Che  Principle  of  Virtual  Work 

The  finite  element  procedure  presented  in  the  preceding  section  is 
formulated  in  the  spatial  domain  and  the  Cauchy  stress  and  the  velocity 
strain  are  employed  to  measure  the  stress  and  strain  staces.  In  this 
section,  the  ALE  finite  element  for  the  equilibrium  equation  is 
formulated  with  a  recourse  to  the  principle  of  virtual  work  in  the 
referential  system.  It  should  be  noted  that  the  linearization 
procedures  in  this  section  are  performed  by  keeping  the  referential 
coordinates  constant  which  is  in  contrast  to  the  usual  race 
formulations  with  the  material  coordinates  held  constant. 

Denoting  a  virtual  variation  of  displacement  by  Che  principle 

of  virtual  work  requires  (Malvern,  1965] 


J  5ui,jTijdRx  "  /  h$ukhkdr*  +  1  5ukbkdRx 

R-x  rx 


(2.25) 


where  the  integration  extends  over  the  spatial  domain.  This  expression 
may  be  transformed  to  the  referential  coordinate  system  as  follows: 
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In  the  above,  Js  is  the  scalar  ratio  of  differencial  areas  drx/dr.(  and 
J  is  the  determinant  of  the  mapping  censor  between  spatial  and 
referential  coordinates  given  by 


Eq.  (2.26)  can  be  recognized  as  the  balance  of  virtual  internal  and 
external  energy  integrated  over  the  referential  coordinates. 
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The  tangent  stiffness  matrix  corresponding  to  Eq.  (2.26)  can  be 
obtained  by  considering  the  rate  form  for  the  virtual  internal  energy 
with  che  referential  configuration  held  constant: 
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As  can  be  seen,  three  rate  quantities  reside  in  Eq .  (2.29).  The  firs 
race  term  can  be  manipulated  by  considering  the  rate  form  for  che 
identity 


3Xm  3xk 
3xj  3Xm  “ 


(2.3 


and  it  cap  be  shown  that 

, 3Xm^  3xk  m  _  3Xq  3 x'< 
x3Xj ‘  3xm  *  3xj  3Xn 
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With  che  substitution  of  mesh  velocity  [Hughes  et  al.,  1973] 


(2.32) 


Eq.  (2.31)  cakes  che  form 
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The  second  race  term  in  Eq.  (2.29)  involves  che  race  of  change  of  the 


Cauchy  stress  tensor  which  has  been  shown  in  Eq.  (2.7) 


"ckTij,k  +  Cijk£v(k,£)  +  Sijkzv[k,*; 


(2.34) 


In  this  equation,  -c,  t  ,  is  the  transport  of  the  stress  histories. 
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The  term  ^  can  be  interpreted  as  the  pure  deformation  part 


for  the  race  of  change  of  che  Cauchy  stress.  The  fourth  order 


generalized  material  response  censor  C^^  is 
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if  Che  Jaumann  stress  race  [Prager,  1961]  is  used ,  whereas 
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if  che  Truesdell  stress  race  is  employed.  The  fourth  order  censor 
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C . . ,  is  given  by 
ijki 
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For  most  of  che  currently  used  material  models,  che  material  response 

c 

censor  C  ^  possesses  both  major  and  minor  symmetries.  Several 

examples  for  che  generalized  material  response  censor  nave  been 

discussed  in  Liu  [1984]  which  will  not  be  repeated  in  che  presenc 

study.  The  term  S  v.  .  in  Eq.  (2.34)  is  the  rocacional  part  for 
ij  *ti  1  ^ ,  2.  J 


the  rate  of  change  of  Cauchy  stress  and  che  fourth  order  tensor  S 


ijkj. 


is  given  by 


Sijki 


2  ^‘ii^jk  'ji^ik  ‘ik^ji  Tjk^ii^ 


(2.38 


The  last  race  quantity  in  Eq.  (2.29)  is  given  by  Liu  [1984] 


J  -  Jv,  , 
k,  X 


Eq .  (2.33)  and  Eq.  (2.39)  can  be  substituted  into  Eq.  (2.29)  to  yield 


y  >  f  3lSul  r  3xm  3^k  .  *_  _  .  3xm  -  |  r.R  (-?  401 

lnC  3Xm  3xk  3Xj  T1^  3xj  Ti''  3Xj  Tl3  <,!C  d  x 


which  can  be  transformed  back  to  the  spatial  domain  as 


wint  "  /  ^Ui,j^Tij  Tij'3k,k  “  Tik^j,kjdE^x 


(2.41) 


In  case  of  v  »  v,  i.e.,  (  )  ■  (’),  Eq.  (2.41)  degenerates  to 


*int 


f  5ui,j^Tij  +  Tijvk,k  ‘  Tikvj,kJdRx 
Rv 


(2.42) 


' x 


which  is  identical  to  the  incremental  expression  for  the  virtual 
internal  work  in  the  Lagrangian  description  [Liu,  1984],  The  rate 
constitutive  equation,  Eq.  (2.34),  may  be  substituted  into  Eq.  (2.41) 
to  yield 


wint  “  /  ^ui, j  ^”ckTij ,k  +  Cijk£v(k,£)  +  Sijk£v[k,i] 
*x 

+  T ij^k.k  "  Tikvj,’JdRx 


When  vi  m  -  ci  is  substituted  into  the  last  two  terms,  Eq.  (2.43) 


can  on.be  written  as 


[* 


winc  *  i  5u-  -D- -t  vi  aR  +  ja.  .1...  v, 

p  4»J  x  ‘d  1 . J  IJKl  <,J 

*-x  *x 


+  f  dui,j(TikCj,k  '  Tijck,k  "  ckTij,k)dR* 


(2.47) 


in  which 


Q  C 

ijkl 


for  Truesdell  rate 


for  Jaumann  rate 


(2.43) 


T 

ijiti 


:  .  o  . , 
J£  ik 


(2.49) 


since  the  Truesdell  stress  rate  has  been  chosen  for  the  Cauchy  stress 
tensor.  Note  that  the  first  two  terms  in  Eq.  (2.47)  are  identical  with 
the  linearized  expression  for  the  virtual  internal  work  in  the 
Lagrangian  description.  The  transport  term  arises  when  cj  *  0  (i.e., 
when  the  ALE  description  is  used)  and  it  is  independent  of  the  stress 
rate  chosen  for  the  Cauchy  stress  tensor. 


Discussions  on  the  Virtual  Internal  Work  in  ALE  Description 

The  domain  of  integration  for  Eq.  (2.47)  C3n  be  transformed  to  the 
material  system  [Liu,  1984]  as 


*  ,  Sou*  3 Xj  3Xi.  ,  3v, 

winc  *  77 —  +  T71  77^  c  •  >  dRX 

Rv  3Xj  1  21  3XP  3Xq  ?jq£  dXl  x 


,  3xt  3  Xu 

+  /  5a.  -7-  ra£  S„  c4  ,  dRy 

d»J  3^p  Pq  ji14 


3  xt  3  x, 

f  5u.  .  y  -  y 1  S  c,  ,  dRy 

J  i,J  3Xp  3Xq  ?q 


fRx5ui  jVij,kdeCWdRx 


(2.50) 


In  this  expression,  the  second  Piola-Kirchhoff  stress  S  is  related  to 

pq 


the  Cauchy  stress  by 


1  3  Xy  3  x 


ij  3x  3X_  3X  pq 

detfrrr)  H  M 


(2.51) 


The  fourth  order  tensor  CT.,  is  the  counterpart  of  the  material 

ijRfc 


response  tensor  for  the  Truesdell  stress  rate  in  the  material  domain 


3x  3  Xy  3X.,  3X^  ^Xj^  r 


CijkZ  3  deC^3X^  3 x_  3x  3xr  3 


ip  3xq  jx.  dxs  pars 


(2.52) 


Ihe  material  velocity  defined  bv 


(2.53) 


v, 

< 


can  be  expressed  by  using  the  chain  rule 


(*) 


C) 


(2.54) 


as 


* 

\ 


3xj.  * 

3X7  Xi 


(2.55) 


The  derivative  of  material  velocity  with  respect  to  the  material 
coordinates  can  be  written  as 


*  it 

3v^ 

3X~  “  lx"  '  Tx“  lx" 

l  l  i  l 


(2.56) 


It  is  noted  that  the  second  derivative  of  spatial  coordinates  with 

* 

respect  to  the  material  coordinates,  -  „ v  X<  is  neglected  at  this 
point  because  only  the  first  derivative  terms  are  required  in  the 
derivation  of  the  tangent  stiffness  matrix. 

The  convective  velocity,  c^  *  vk  ”  *  can  expressed 

in  terms  of  the  rate  of  material  and  spatial  coordinates  as 


in  which  Che  expressions  for  che  material  and  mesh  velocities  given  by 
Eq.  (2.55)  and  Eq.  (2.32),  respectively,  have  been  used.  The  spatial 
derivative  of  che  convective  velocity  can  be  shown  as 


3  xu  * 

c,  .  -  -  ~  X.  . 
k,j  3Xt  1,2 


(2.58 


Eqs.  (2.56-58)  can  be  substituted  into  the  race  of  change  of  virtual 
internal  work,  Eq.  (2.50),  to  yield 


f  c  rik  3xk  - 

iat  J  aXj  sijl-3Xi  '  3Xn  sxj^x 


Rx 


★  * 
36uu  3 Xu  3x  ,3x,.  3x^  3X 

+  /Rx  ixj"  3X7  3X7  CijPq  3X7  ‘  3*7 


36uk  3xk  3Xn 
[RX  3^  *X1  U  X 


3<5uk  3xv  3Xn  ^ 
JRx  3Xj  3Xi  ^ij  3Xn  dRX 


3x. 


(2.59 


It  is  noted  that  che  first  integral  in  Eq.  (2.59)  represents  the 
initial  stress  effect;  che  second  incegral  is  the  material  response 
part;  the  third  and  fourth  integrals  arise  because  of  the  use  of  ALE 


description  (c^  ?*  0);  and  the  last  incegral  represents  the  transport  of 
the  stress  histories  in  the  ALE  description. 

An  alternate  derivation  of  Eq.  (2.59)  is  presented  in  Haber  and 
Abel  [1983].  This  formulation  is  subsequently  improved  by  Haber  [1984] 
in  which  the  so-called  "Mixed  Eulerian-Lagrangian"  description  is 
inspired  by  the  work  of  Hughes  et  al.  [1978].  In  this  paper,  the 
second  Piola-Kirchhof f  stress,  the  Green  strain,  and  the  virtual  work 
expression  are  linearized  to  obtain  the  tangent  stiffness  matrix. 
Because  of  the  choice  of  this  conjugate  pair  of  stress  and  strain,  and 
the  linearizaton  procedure,  tedious  algebra  is  involved  and  the 
expression  for  the  tangent  stiffness  matrix  is  only  valid  for  linear 
elastic  isotropic  material.  For  comparison  purposes,  the  rate 
counterpart  of  the  Haber's  formulaton  is  re-derived  in  Appendix  B.  The 
differences  between  the  present  and  the  Haber's  formulations  are 
summarized  in  Table  2.1.  As  can  be  seen,  the  transport  of  stress 
histories  is  not  included  in  Haber's  formulation.  This  deficiency 
results  from  the  choice  of  the  second  Piola-Sirchhof f  stress  and  the 
Green  3train,  and  the  linearization  procedure  for  the  constitutive 
equation.  The  transport  of  stress  histories  is  very  important  for 
highly  convective  calculations  such  as  the  wave  propagation  problems 
presented  in  the  next  chapter. 


Table  2.1.  Comparisons  between  Che  Present  and  the  Haber's 
Formulations 


Present  Formulation 


Haber  [1984] 


stress  t , .  (Cauchy  stress) 

ij 


S  (second  Piola-Kirchhof f 
3  stress) 


strain  v  .  (velocity  strain) 

1  >  J 


E,  .  (Green  strain) 

ij 


constitutive  t  -  -V  *  cIju\k,t)  sij  '  'ijkAl 
equation  J  J 


*  SijklV[k,2] 


where  C.  ,  ■  constant 

ljk  l 


expressions  36uj,  A  ^ 

for  ;1K  fRx  3Xi  '  sx«  3Xi]  % 


,  3<5uk  3xk  3xr  _s  <•  3xr  _  3xr  3‘tn,  . 

+  Ry  3Xj  3Xi  3Xq  CiJP(’'3Xp  3Xn  3Xp;  X 


same,  except 

C® .  is 

ijpq 


replaced  by 


This  term  is  omitted 


computer  Implementation  of  the  ALE  Tanzenc  Stiffness  Matrices 


Introducing  the  interpolating  functions 


5ui  *  ^a^uia  *  a  ■  1,  number  of  nodes  per  element  (MEN)  (2.60) 


,  b  -  1,  MEN 


(2.61) 


the  matrix  equation  associated  with  the  principle  of  virtual  work  reads 


Ktan  ~  „  fext  .  fint 


(2.62) 


The  element  tangent  stiffness  matrix  corresponding  to  the  race  form  of 
virtual  internal  work,  Eq.  (2.47),  is  composed  of  three  matrices 


Ktan  „  kD  +  kG  +  KALE 


(2.63) 


In  computer  implementation,  it  is  inconvenient  to  deal  with  the 
indices  of  a  fourth  order  tensor  such  as  Hence,  the  procedure 

given  in  Liu  and  Ma  [1982]  will  be  followed  to  develop  the  element 
tangent  stiffness  matrices.  Interested  readers  are  suggested  to  refer 
to  this  work  for  the  expressions  of  and  which  will  not  be 
repeated  here.  The  matrix  can  be  shown  as: 


K  ,  -  f  Bi  3  dR  -  f  3iT,;3vdR 

— ab  ‘  ^  ~a~I~b  x  1  ^  ~a~2~b  x 


where  Che  9x3  matrix  B,  is  given  bv 

~b 


Vl 

0 

0 

^b,2 

^b,i 

0 

0 

*b,2 

0 

0 

0 

*b,3 

0 

Bh  " 

~b 

^b,3 

*b,3 

0 

*b.2 

V 

*b,2 

‘*b,l 

0 

0 

*b,3 

"*b,2 

~*b,3 

0 

*b,l 

The  3x3  matrix  is 


~b 


0 

9, 


0 

0 

9> 


'he  9x9  unsymmecric  matrix  T,  is 


The  9x3  T-  matrix  involves  the  spatial  derivatives  of  the  Cauchy  stress 
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Tl,l 

T  1,2 

T  1.3 

t2,1 

t2,2 

T  2 , 3 

T3 , 1 

t3,2 

T  3 , 3 

T4,l 

T4,2 

T*,3 

T5',l 

T5.2 

T5,3 

T6,l 

T6,2 

T6,3 

0 

0 

0 

0 

0 

0 

0 

0 

0 

(2.63) 


These  ALE  tangent  stiffness  matrices  may  be  incorporated  into  any 
existing  Lagrangian  finite  element  code  without  changing  the  main  body 
of  Che  program,  nevertheless,  the  amount  of  computations  is  scill 
substantial.  For  this  reason,  the  ALE  tangent  stiffness  matrices  are 
further  rearranged  as  the  product  of  several  3x3  matrices  as  follows. 

The  ALE  tangent  stiffness  matrix  corresponding  to  Eq .  (2. -*7)  can 


be  arranged  in  three  parts  as 


The  first  matrix  can  be  shown  as 


S*  ■  ! 


(2.70) 


where  che  components  of  Che  3x3  t  matrix  are 


T1 

T2 

t6 

I  -  • 

t2 

T  3 

t5 

(2.71) 

t6 

T  5 

T4 

components  of 

the 

3x3  £b 

matrix  are  given  by 

*b,l 

*b,l 

*b,l 

ib  "  >bj2 

♦b.2 

*b,2 

(2.72) 

*b,3 

*b,3 

*b,3 

and  che  components  of  Che  3x3  $  ^  matrix  are 


>a,l  °  ° 

°  S  a ,  2  ° 


(2.73) 


0 


3 


The  second  matrix  is  given  by 


(2)  - 
kV  -  /  .  t*  i’dR 

ab  J  n  — ‘-a^b  x 


and  Che  lasc  matrix  K^)  is 


£{Vm  /  T’i.dR 
~ab  J  _  ~a~b  x 

‘Sc 


where  the  components  of  the  3x3  matrix  are 


t’  “  $ . 


Lj,l  '  1  j  .2  *  Ij  .3 


~a  a,  j  2j ,  1  2j  ,2  ‘2j,3 

T3j,l  T 3j  ,2  T 3j , 3 


The  components  of  the  matrix  4,  have  already  been  given  by  Ec 


(2.75) 


.  (2.66). 
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CHAPTER  3 

.NUMERICAL  PROCEDURES  FOR  ALE  EQUATIONS 

3.1  Overview 

The  ALE  formulation  for  conservation  laws,  constitutive  equations, 
and  equation  of  stace  have  been  considered  and  derivations  from  strong 
forms  to  matrix  equations  have  also  been  given  in  the  previous  chapter. 
However,  the  overall  effectiveness  of  che  present  effort  depends  to  a 
large  degree  on  the  numerical  procedures  used  for  the  system  equations. 
It  is  the  purpose  of  this  chapter  to  present  some  numerical  methods  for 
the  computations  of  the  ALE  equations  which,  in  general,  are  impossible 
to  integrate  analytically. 

Among  the  vast  amount  of  numerical  methods,  subjects  relevant  to 
ALE  computations  include  the  evaluations  of  generalized  mass  matrices, 
internal  and  external  force  vectors,  nonlinear  convection  effects,  and 
new  stress  and  strain  states  in  addition  to  the  construction  of  time¬ 
stepping  algorithms.  Since  discussions  or  the  mass  matrix  and  force 
vectors  have  been  well-documented  in  the  literature  [Zienkiewicz , 

1977],  these  treatments  will  not  be  included  in  this  investigation. 

For  the  on-growing  area  of  nonlinear  convective  computations,  a  review 
of  various  approaches  such  as  upwind  methods,  Pecrov-Galerkin ,  and 
Tav lor-Gaierkin.  formulations  is  provided  in  the  next  section.  The 
computations  of  scress  and  strain  states  for  path-dependent  materials 
3re  then  presented  in  Section  3.3,  which  includes  the  discussion  of  the 
stress-velocity  product  and  a  collocation  weighted  residual  formulation 


for  Che  constitutive  equations.  In  Section  3.4,  the  nonlinear 
computational  procedures  for  a  two  quadrature  point  element  is 
presented.  An  explicit  algorithm  for  the  integration  of  the  ALE 
equations  is  shown  in  Section  3.5  and  the  application  of  the  present 
formulation  to  a  number  of  problems  is  provided  in  Section  3.6. 


3.2  Finite  Element  Methods  for  Non-Self-Adjoint  Equations 


The  convective  effect,  which  arises  when  the  (Quasi-)Eulerian 
description  is  employed,  has  been  one  of  the  difficult  topics  in  the 
development  of  numerical  methods.  The  usual  ( Bubnov-)Galerkin  method 
leads  to  a  non-symmetric  convective  operator,  and  spurious  spatial 
oscillations  are  exhibited  at  moderate  to  high  Peclet  numbers.  Though 
highly  refined  meshes  can  remove  these  oscillations,  the  advantages  of 
the  finite  element  method  may  be  diminished. 

Several  remedies  have  been  devised  to  overcome  this  difficulty. 

In  the  finite  difference  category,  the  use  of  upwind  differencing  on 
the  convective  term  is  proposed  in  Richtmyer  and  Morton  [1967]  and 
later  detailed  in  Spalding  [1972].  A  similar  idea  is  introduced  into 
the  finite  element  method  in  Christie  et  al.  [1976]  in  which  the  weight 
functions  are  skewed  to  achieve  the  upwind  effect.  Relevant  early 
articles  in  these  upwind-type  finite  elements  include  Heinrich  et  al. 
[1977],  Hughes  [1978]  and  3elytschko  and  Kennedy  [1978],  among  others. 
Of  particular  interest  in  the  present  ALE  calculations  is  the  work, 
presented  in  Brooks  and  Hughes  [1982].  In  this  paper,  the  artificial 
diffusion  operator  is  constructed  in  a  tensorial  form  so  as  to  act  only 
in  the  flow  direction  (streamline  upwind).  The  free  parameter  In  this 
method  is  the  amount  of  diffusion  selected  to  maximize  the  solution 
accuracy.  Detailed  discussion  of  this  streamline  upwind  method  is 


beyond  the  scope  of  the  present  investigation;  nevertheless, 


derivations  of  this  method  for  the  Navier-Stokes  equation  in  the  ALE 


description  can  be  found  in  Liu  [1980],  and  this  method  will  be 
employed  to  handle  the  convective  effect  at  the  present  stage. 

In  addition  to  the  above  school  of  upwind  techniques,  the 
Petrov-Galerkin  method  proposed  in  Dendy  [1974]  consists  of  choosing 
different  classes  of  functions  for  the  weight  and  trial  functions. 
Because  of  the  potential  of  Petrov-Galerkin  methods  for  flow  problems, 
a  large  amount  of  literature  has  been  accumulated.  Among  various 
techniques,  the  work  presented  in  Hughes  and  Tezduyar  [1984]  and  Morton 
et  al.  [1980]  appear  adapted  to  deal  with  the  present  ALE  formulation. 
It  is  also  worthwhile  to  point  out  that  the  Taylor-Galerkin 
formulations  given  in  Donea  [1984]  and  Lohner  at  al.  [1984]  appear  to 
be  potential  candidates  to  handle  the  convective  effect  encountered  in 


the  ALE  computations 


3«3  Stress  Update  Procedures 

The  stress  state  of  a  path-dependent  material  depends  on  the 
stress  history  of  the  material  point.  A  stress-history  can  be  readily 
treated  in  a  Lagrangian  description  because  elements  contain  the  same 
material  points  regardless  of  the  deformation  of  the  continuum; 
similarly,  quadrature  points  at  which  stresses  are  computed  in 
Lagrangian  elements  coincide  with  material  points  throughout  the 
deformation.  On  the  other  hand,  in  an  ALE  description,  a  mesh  point 
does  not  necessarily  coincide  with  a  material  point  so  that  the  stress 
history  needs  to  be  convected  by  the  relative  velocity  c,  as  indicated 

A# 

in  Eq.  (2.7).  Note  that  the  spatial  derivatives  of  the  deviatoric 
stress  are'  involved  in  the  convection  term. 

When  C-1  functions  are  used  to  interpolate  the  element  stresses, 
the  ambiguity  of  the  stress  derivatives  at  the  element  interface 
renders  the  calculation  of  the  spatial  derivatives  of  stress  a 
difficult  task.  As  mentioned  in  REMARK  2.1.3,  this  is  remedied  by 
replacing  Eq.  (2.7)  by  a  set  of  coupled  equations,  Eqs.  (2.11-12), 
and  the  corresponding  matrix  equations  have  been  given  in  Eqs. 
(2.22-23).  It  should  be  noted  that  all  the  path-dependent  material 
properties,  such  as  yield  strains,  effective  plastic  strains,  yield 
stresses  and  back  stresses,  should  be  convected  via  this  procedure  with 
s  replaced  by  each  of  these  properties  in  turn,  and  with  z 
appropriately  modified. 

In  this  section,  numerical  methods  for  the  definition  of  the 
stress-velocity  product  and  computation  of  the  incremental  stresses 


are  presented.  Although  the  following  development  is  discussed  in  a 
two-dimensional  setting,  the  extension  to  the  three-dimensional  case  is 
straightforward. 


Construction  of  the  Stress-Velocity  Product 

In  a  nonlinear  displacement  finite  element  formulation,  the 
velocities  are  stored  at  nodes  while  the  stress  histories,  back 
stresses  and  yield  radii  are  available  only  at  quadrature  points.  In 
order  to  establish  the  nodal  values  for  the  stress-velocity  product,  a 
weak  formulation  is  a  logical  necessity.  In  addition,  based  on  the 
one-dimensional  study  presented  later  in  this  chapter,  in  which  the 
upwind  procedure  is  used  to  define  this  intermediate  variable,  the 
artificial  viscosity  technique  (streamline  upwind)  [Brooks  and  Hughes, 
1982]  is  considered  here  as  a  generalization  of  this  upwind  procedure 
to  multi-dimensional  cases.  For  the  sake  of  clarity,  the  free  indices 
i  and  j  denoting  the  component  of  stress  tensor  will  be  dropped 
hereafter. 

The  relation  for  the  stress-velocity  product  given  in  Eq.  (2.12) 
is  modified  to  accommodate  the  artifical  viscosity  censor  A^ 

7k  -  sck  -  Akm,m  (3>1) 

The  ingredients  of  the  artificial  viscosity  tensor  consist  of  a 


Censorial  coefficient  multiplied  by  Che  stress: 


where  Che  censorial  coefficienc  is  conscrucced  co  acc  only  in  che  flow 
direccion  (screamline  upwind  effecc)  [Brooks  and  Hughes,  1982] 


ukm  *  uCkca/cncn  (3.3) 

and  che  scalar  y  is  given  by 

NSD 

U  -  C  c^c^/NSD  (3.4) 

i-1 

Here  h^  is  che  elemenC  lengCh  in  Che  i-direcCion,  NSD  designaces  che 
number  of  space  dimensions,  and  a*  is  che  arcificial  viscosicy 
parameCer  discussed  in  Appendix  B  given  by 

0^  *  ~  y  ,  for  c^  \  0  (3.5) 

The  weak  form  corresponding  Co  Eq.  (3.1)  can  be  obcained  by 
mulciplying  by  che  cesc  funccions  for  che  scress-velocicy  produce  and 
incegracing  over  che  spacial  domain  R^. 


/  ‘VA  ■  I  sV\dK*  - 1  ‘Vta.,."* 

^x  ^ 


This  equacion  may  be  wriccen  as 


I  5VkdR* 

*x 


{  5yksckdR* 
^x 


5yk,mAkmdR^ 

‘Sc 


(3.7) 


by  applying  the  divergence  theorem  and  by  assuming  no  traction  asso¬ 
ciated  with  the  artificial  viscosity  on  the  boundary.  The  expression 
for  Aka,  Eqs.  (3.2-3),  can  be  substituted  into  this  equation  to  yield 


/  +  6K)3CkdRx 

*x 


(3.8) 


where 


Syk  -  'Syk>ai;icin/cncT1  (3.9) 

can  be  viewed  as  a  modification  of  the  Galerkin  finite  element  method 
because  of  the  transport  nature  of  stress-velocity  product. 

The  shape  functions  for  the  stress-velocity  product  can  be  chosen 
to  be  the  standard  C°  functions.  The  number  and  position  of  numerical 
integration  points  for  Eqs.  (3.8-9)  should  be  selected  to  be  the 
quadrature  points,  since  the  stress  histories  in  Eq.  (3.3)  are  only 
available  at  these  points. 

Remark  3.3.1:  The  determination  of  the  artificial  viscosity  parameter 
is  accomplished  with  a  recourse  to  the  analytic  solutions  for  a 
one-dimensional  constitutive  equation.  The  analytic  solutions  are 
obtained  by  employing  the  Laplace  Transform  technique  which  is  detailed 


in  Appendix  C 


Multiple  Stress  Collocation  Point  Formulation 

Following  the  procedures  given  in  the  previous  sub-section,  the 
stress-velocity  product  can  be  defined  at  each  nodal  point  and  it  can 
be  substituted  into  the  constitutive  equation  as  follows  to  calculate 
the  rate  of  change  of  stresses. 

* 

Mss  »  z  -  G^y  +  D  s  (3.10 

Definitions  for  the  above  matrices  and  vectors  are  given  in 
Appendix  A.  Note  that  the  interpolation  functions  for  stresses  need 
to  be  integrated  over  the  spatial  element  domain  in  these  definitions. 
However,  the  present  displacement  formulation  carries  stresses  only  at 
quadrature  points  in  contrast  to  the  Hellinger-Reissner  and  Hu-Washizu 
[Washizu,  1975]  finite  elements  in  which  the  stress  interpolants  are 
assumed . 

Furthermore,  the  task  to  select  the  number  of  quadrature  points 
for  the  displacement  finite  element  poses  another  important  issue.  Foi 
example,  the  locking  phenomenon  for  fully  integrated  elements  arises 
when  the  material  becomes  incompressible.  While  selective  reduced 
integration  can  overcome  this  difficulty,  it  is  just  as  costly  as  full 
quadrature.  To  alleviate  this  computational  hurdle,  the  use  of  one 
point  quadrature  combined  with  hourglass  control  is  developed  in 
Belvtschko  et  al .  [1984],  In  addition,  a  nonlinear  two-quadrature 


point  element  presented  in  the  next  section  appears  to  be  another 
candidate  for  large  scale  computations  because  it  exhibits  nearly  the 


^  ^  un  d*  *.^»  . 


same  accuracy  as  che  selective  reduced  integration  element  while  with 
only  one-third  of  the  cost.  The  elements  mentioned  here,  as  well  as 
others  in  the  family  of  displacement  elements,  can  be  readily  adopted 
in  the  ALE  computations.  It  is  obvious  chat  a  stress  transport 
procedure  suitable  for  any  number  of  quadrature  points  is  needed. 

Inspired  by  the  equivalence  proof  for  the  mixed  and  displacement 
elements  in  Belytschko  et  al.  [1985],  che  displacement  element  is 
divided  into  M  subdomains  where  M  denotes  the  number  of  quadrature 
points.  Each  subdomain  is  designated  by  Rj  (I  *  1,  M) ,  which  contains 
the  quadrature  point  x^,  and  no  two  subdomains  overlap.  Associated 
with  Rj,  a  stress  interpolating  function  is  assigned  and  its  value 
is  prescribed  only  at  quadrature  point  x  *  x_  to  be  unity,  or 


n^i5  *  1 


The  test  function  in  Rj  is  chosen  to  be  the  Dirac  delta  function 


(3.11) 


♦!  *  5(x-xi) 


(3.12) 


Substitutions  of  these  functions  into  the  constitutive  equation  given 
by  Eq.  (2.16)  represents  a  mathematical  requirement  that  che  residual 
of  the  weak  form  vanishes  at  each  collocative  quadrature  point. 

Because  the  collocation  point  is  located  right  at  the  quadrature  point, 
the  algebraic  equations  resulting  from  Eq.  (3.10)  are  dependent  only  on 
the  information  (stress  history)  associated  with  its  host  points. 

Since  the  Dirac  delta  function  has  che  important  property  chat 


(3.13 


/  F(x)<5(x-Xj)dRx  -  F(xj) 

*x 

All  of  Che  matrices  and  vectors  in  Eq.  (3.10)  can  be  easily  worked  out 
without  numerical  integration  and  they  are  given  in  Appendix  D. 

Remark.  3.3.2:  In  Appendix  E,  the  stress  update  procedure  is  given  for 
a  uniform  one-dimensional  mesh;  the  resulting  stress  convective  cerms 
bear  a  strong  similarity  to  the  donor-cell  differencing  [Roache,  1972] 
On  the  other  hand,  a  simple  averaging  (or  the  central  differencing)  is 
obtained  if  the  upwind  effect  is  not  applied  to  the  convective  terms. 


3.4  An  Efficient  Nonlinear  Element 

In  Liu  et  al.  [1985],  a  "unified”  stabilization  procedure  for  the 
Laplace  equation,  the  equation  for  continuum  and  the  equation  for 
plates  and  shells  in  two  and  three  dimensions  is  presented.  The  major 
achievements  are  (l)  control  of  the  spurious  singular  modes  for 
underintegrated  finite  elements;  (2)  enhancement  of  the  computational 
efficiency  without  sacrificing  the  accuracy;  and  (3)  alleviation  of  the 
locking  phenomenon  for  continuum  elements  when  the  material  becomes 
incompressible.  Accurate  solutions  are  achieved  for  linear  isotropic 
(incompressible)  materials.  However,  when  the  above  IPS  element  is 
applied  to  plastic  materials,  the  solution  accuracy  deteriorates  and 
Che  proposed  element  tends  to  be  too  stiff  when  compared  to  the 
selective  reduced  integration  (SRI)  element. 

A  possible  reason  for  this  shortcoming  in  nonlinear  calculations  is 
due  to  too  few  (only  one)  stress  sampling  points  in  each  element.  The 
entire  element  must  be  in  either  a  purely  elastic  or  a  purely  plastic 
state  determined  by  evaluating  the  constitutive  equations  at  the 
element  center.  Consequently,  the  onset  of  plastic  fronts  and  the 
plastic  yielding  effects  cannot  be  accurately  accounted  for.  On  the 
other  hand,  the  SRI  element  permits  different  stress  states  in  elements 
because  the  constitutive  equations  are  calculated  at  four  (for  a 
two-dimensional  bilinear  interpolation  of  displacement)  quadrature 
points.  Nevertheless,  in  view  of  the  large  additional  cost  in 
compuiing  Che  spatial  derivatives  of  shape  functions  and  in  evaluating 


the  constitutive  equations,  there  seems  to  be  little  benefit  in  using 
the  SRI  element. 


Based  on  these  considerations,  research  effort  has  been  directed 
toward  the  development  of  an  accurate  two  stress-sampling  point  element 
(IPS2),  which  retains  the  efficiency  of  the  IPS  element.  The  basic 
idea  of  the  present  approach  is  to  adopt  the  approximation  of  the 
gradient  operator  matrix  in  Liu  et  al.  [1984]  and  to  evaluate  the 
constitutive  equation  at  two  presumed  quadrature  points.  For 
illustrative  purposes,  details  of  the  two-dimensional  plane-strain  case 
are  presented  here. 

The  gradient  operator  matrix  given  in  Liu  et  al.  [1984]  is 


5a 


Ba(0) 

Al 


♦  Bf>>„ 


(3.14) 


where  the  subscript  ''a”  ranges  from  1  to  the  number  of  nodes  per 
element.  Definitions  of  B,(0),  B^^CO)  and  Bde^(Q)  can  be  found  in 

**  a*®  » n  >n# 

Appendix  F.  The  incremental  strain  is 

Ae($.n)  -  (3,(0)  +  Bj*I(0)«  +  BaSn(0),l^da  (3.15) 

AM  AM  <V°  I  AM  Al“  9*1  A»  Al“ 

where  Ad,  is  the  incremental  displacement  vector.  The  incremental  spin 
censor  can  similarily  be  approximated  by 


AW(e,n) 


£a,n(£)nlida 


(3.16) 


where  WA(0),  W  ^(0)  and  V  (0)  are  given  in  Appendix  G.  The 
incremental  strain  and  spin  tensors  are  calculated  at  two  quadrature 
points , 


✓3  ’  '  /3; 


1_  -  i_v 

✓3  *  +  /3 


(3.17a) 


(3.17b) 


and  the  stress  history,  xn^c,j),  and  the  back  stresses,  an(?i),  at  time 
step  n  are  rotated  by  Q(£j2>)  Co  account  for  the  rotational  effect 
[Hughes  and  Winget,  1980]. 


Xn+l^az^  ^2  In  2i^i 


(3.18a) 


(3.13b) 


where  Tn>i(^£)  and  are  the  intermediate  Cauchy  stresses  and 

back  stresses,  respectively.  The  expression  for  the  orthogonal  matrix 


0  -  I  +  (I  -  ^  iW)_1AW 

~  •—  L  — 


(3.19) 


where  I  is  the  identity  matrix.  The  incremental  change  of  the  stress 
state  due  to  material  deformation  at  the  quadrature  point  r  ,  can  be 


handled  by  che  plasticity  integration  procedure.  An  example  is  the 
radial  return  method  given  in  Krieg  and  Key  [1976]. 


*!<£*>  *£t(a£*  £n+l*  £n+l »  kn,-*‘>  (3‘20) 

M  £j(a£>  £n+l  *  £n+l  *  ’S,*”)  C3‘21) 

£l>  *  fk(A£’  fn+1*  £n+l »  kn,“-»)  (3,22) 


where  Ak  is  the  incremental  change  of  the  yield  radius;  f  ,  ,  and 

fj,  designate  che  plasticity  integration  procedure  for  at,  Ao,  and  Ak, 
respectively.  The  state  of  stresses  at  step  ni-l  is  then  updaced  by 

Xn+l^Z3  *  In+l^  +  (3,2 

En+l^lz5  ”  £n+l(£z)  +  A£^£z)  (3'2 

^n+l <£i>  "  kn<£z}  +  Ak(£i>  C3’2 

The  element  internal  force  vector  is 


1  -  /  iTIn.l  (3-26) 

The  approach  used  in  Eq.  (3.14)  can  be  applied  again  to  approximate  the 
gradient  operator  matrix  B  in  Eq.  (3.26).  Numerical  integration,  with 
weights  equal  to  2  and  with  the  Jacobian  assumed  to  be  a  constant. 


e>: 
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(element  area)/4,  is  employed  to  integrate  the  element  internal  force 
vector 


7  +5a "<£>«, 
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(3.2: 


and  A  represents  the  area  of  the  element. 

The  element  internal  force  vector  can  be  further  arranged  as 


£  “  £l  +  £stab 


(3.2! 


where  the  internal  force  resulting  from  the  one  point  quadrature  is 


<Tll>3l  +  ^12^2 

<t22>£2  +  <t12>^1 


(3.2< 


and  the  contribution  from  the  stabilization  procedure  is 


£stab  ”  e 
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In  these  expressions,  the  following  definitions  have 
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Table  3.1 


Comparison  of  Computations  and  Storage  for 
SRI ,  IPS  and  IPS2  Elements 


2-D 

Example 

Spatial  Derivatives 
of  Shape  Functions 

Evaluation  of 
Constitutive 
Equation 

Storage 

Requirements 

SRI 

4  times 

4  times 

(3C.S.+38.S.+1Y.R.) 

*4  «  28  words 

IPS 

1  time 

1  time 

9C.S.+3B.S.+1Y.R. 

-  13  words 

IPS2 

1  time 

2  times 

(3C.S.+3B.S.+1Y.R.) 

*2  *  14  words 

C.S.:  Cauchy  Stress;  B.S.:  Back  Stress;  Y.R.:  Yield  Radius 


Foe  simplicity,  the  coupled  equations  will  be  integrated  by  an 


explicit  scheme.  Lumped  mass  matrices  are  used  to  enhance  the 
computational  efficiency.  If  (  )n  and  (  )n+1  denote  the  matrices  at 
times  tn  ■  nAt  and  tn+^  ■  (n+l)At,  respectively,  where  At  is  the  time 
increment,  the  explicit  predictor-corrector  mechod  [Hughes  and  Liu, 
1978]  gives 
The  mass  equation: 


+  &+1) 


£n+l  "  £n  +  d-oJAt^ 


an+1  *  £n+l  +  aAtg^i 


The  momentum  equations: 


N  v  .  . ) 
~tv~n+l 


Za  +  (l-^)AcXn 


2n+l  “  Xn+1  +  ia+1 


£q.  (3.37)  needs  to  be  used  in  conjunction  with 


in+1  "  £a  +  *cXn  +  %  "  S)At^vn 


ln+1  *  i+l  +  Sit2  * 


(3.4. 


to  calculate  the  finC. 

~a 


The  energy  equation: 


£n+l  “  ^2n^  ^£n+l  In  Jln£n+P 


(3.4: 


£n+l  "  «n  +  <1  -  C)At«, 


(3.4: 


Sn+1  "  Srt-1  +  Wt  £n+l 


(3.44 


The  constitutive  equation: 


Zn+1  “  (*2n^  12n£n+l 


(3.45 


£n+l  “  ^^n^  ^ (£n+l  ~  £n£a+l  +  £n®a+l^ 


(3.46 
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(3.47 


3-j.i  +  <At  3 


(3.48 


where  a,  8,  y,  ;  and  <  are  the  computational  parameters.  For  explicit 
calculations,  the  following  constraints  on  the  parameters  are  used: 
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a-0;3“0,Y>y;;-0;«:,B0  (3.49) 

The  flowchart  of  the  computational  procedure  for  the  class  of  pressure- 
insensitive  materials  Is  as  follows: 

(1)  Initialization.  Set  n  -  0,  input  initial  conditions. 

(2)  Time  stepping  loop,  t  e  fO,t„,a.,] » 

(3)  Integrate  the  mesh  velocity  to  obtain  the  mesh  displacement  and 
spatial  coordinates. 

(4)  Calculate  incremental  hydrostatic  pressure  by  Integrating  Eqs. 
(3.45-46)  with  s  and  z  replaced  by  p  and  u,  respectively. 

a.  The  rate  of  pressure  due  to  convection. 

b.  The  rate  of  pressure  due  to  deformation. 

(5)  Calculate  incremental  deviatoric  stresses,  yield  scresses  and  back 
stresses  by  integrating  Eqs.  (3.45-46). 

a.  The  rate  of  stresses  due  to  convection. 

b.  The  rate  of  stresses  due  to  rotation. 

c.  The  race  of  scresses  due  to  deformation. 

(6)  Compute  the  internal  force  vector. 

(7)  Compute  the  acceleration  by  the  equations  of  mocion,  Eq.  (3.37). 

(8)  Compute  the  density  by  equation  of  mass  conservation,  Eq.  (3.34). 

(9)  Compute  the  internal  energy  by  the  equation  of  energy 
conservation,  Eq.  (3.42). 

(10)  Integrate  acceleration  to  obtain  velocity. 

(11)  If  (a+l)At  >  t_ax,  stop;  otherwise  replace  n  by  ari  and  go  to  (2). 


p^rrwrriTTTTTwiTTiiTTTT»^rrrrf mibi  Tnir  . . . .  mi  —  M  . 
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3«6  Sumer leal  Examples 

3.6.1  One-Dimensional  Wave  Propagation  Problem 

An  elastic-plastic  wave  propagacion  problem  Is  used  Co  assess  chls 
ALE  approach  in  conjunction  with  Che  simple  averaging  and  Che  proposed 
stress  update  method.  The  problem  statement,  given  in  Fig.  3.1, 
represents  a  one-dimensional,  infinitely  long,  elastic-plastic 
hardening  rod.  Constant  density  and  isothermal  conditions  are  assumed 
to  simplify  the  problem.  It  should  be  noted  that  this  elastic-plastic 
wave  propagation  problem  does  not  require  an  ALE  mesh  and  the  problem 
was  selected  because  it  provides  a  severe  test  of  the  stress  update 
procedure  and  because  of  the  availability  of  an  analytic  solution.  The 
problem  is  solved  using  400  elements  which  are  uniformly  spaced  with  a 
mesh  size  of  0.1.  The  mesh  is  arranged  so  that  no  reflected  wave  will 
occur  during  the  time  interval  under  consideration.  Material 
properties  and  computational  parameters  are  also  depicted  in  Fig.  3.1. 
Four  stages  are  involved  in  this  problem: 

(1)  t  e  [0,  til,  the  mesh  is  fixed.  A  square  wave  is  generated  at  the 
origin. 

(2)  t  e  [ti»  t£ 1 t  the  mesh  is  fixed  and  the  wave  travels  along  the 
bar. 

(3)  t  c  (t2>  C3 ] ,  two  cases  are  studied. 

CASE  A:  The  mesh  is  moved  uniformly  to  the  left  hand  side  with  a 
constant  speed  -v  . 

CASE  B:  Same  as  CASE  A  except  the  mesh  is  moved  to  the  right. 
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1  D  ELASTIC  -PLASTIC)  WAVE  PROPAGATION 


p= 1  E=l(r  E/Et=3  ry= 75  t  =-100 
As=Ax=0.1  v*=0.25VE/p  /?= 0.0  7=0.6 

1.  t€[0,  til  mesh  fised,  wave  generated 


2.  t€[t«,  tg]  mesh  fixed,  wave  travelling 


3A.  t€[t2,  tg]  CASE  A  :  move  mesh  with  v  =  -v 

tj=45  t2=80  1^=160  (10~3) 


x=0 


A  A, 

3B.  t€[t2,  tj]  CASE  B  :  move  mesh  with  v  =  +v 

tj=45  t2=240  tg=320  (10-3) 


x=G 

4.  t  =  tj  report  stress  vs.  spatial  coordinate 


Fig.  1 


(4)  c  *  C3,  che  stress  is  reported  as  a  function  of  spatial 

coordinates  in  Fig.  3.2  and  Fig.  3.3  for  CASE  A  and  CASE  B, 
respectively. 

For  both  cases,  the  momentum  and  stress  transport  are  taken  into 
account  by  employing  the  full  upwind  method  for  elastic  and 
elastic-plastic  materials.  The  results  are  compared  to 

(1)  the  simple  averaging  runs,  in  which  Che  momentum  transport  is 
handled  by  the  full  upwind  method  and  the  stress  transport  is 
computed  by  che  simple  averaging  method,  and 

(2)  fixed  mesh  runs,  in  which  che  finite  element  mesh  is  fixed  in 
space  and  the  results  are  pretty  close  to  the  analytic  solutions. 
The  relative  velocity,  c  ■  v  -  v,  in  CASE  A  (  £  <  0)  is  greater 

than  CASE  B  (v  >  0);  therefore,  the  transport  effect  of  the  former  is 
more  significant.  These  phenomena  have  been  studied  in  Figs.  3.2-3  by 
varying  the  time  step  reported  in  Table  3.2.  The  wave  arrival  time  for 
both  the  proposed  method  and  the  simple  averaging  method  agree  well 
with  the  fixed  mesh  runs.  However,  the  averaging  method  causes  severe 
unrealistic  spatial  oscillations  in  CASE  A  because  of  the  significant 
transport  effect.  The  method  proposed  here  eliminates  these 
oscillations  completely.  Based  on  these  studies,  it  is  found  chat  the 
transport  of  stresses  as  well  as  yield  stresses  (and  back  stresses  for 
kinematic  hardening)  plays  an  important  role  in  the  ALE  computations 
for  path-dependent  materials,  and  che  proposed  update  procedure  is 
quite  accurate  and  effective. 


Table  3.2.  The  seep  sizes  and  numbers  of  time  seeps  for 
elastic-plascic  wave  propagation  example 


• 

*  Time  Step 

At 

Number 

of  Time  Steps 

j 

Cr 

CASE  A 

CASE  B 

56 
0.072 


C  *  Ax/(/E/p  +  |  c  |  ) 


.5 

400 

80 

.7 

286 

57 

.9 

222 

44 

3.6.2  Two-Dimensional  Elastic-Plastic  Wave  Propagation 


In  addition  to  the  preceding  problem,  a  two-dimensional 
plane-strain  elastic-plastic  wave  propagation  problem  is  considered 
here  as  an  another  test  for  both  the  IPS2  element  and  the  proposed 
multi-stress-point  transport  procedure.  A  100x50  mesh  is  used  to  model 
a  spatial  domain  of  size  10x5.  The  radial  return  procedure  given  in 
Krieg  and  Key  [1976]  is  used  and  isotropic  hardening  is  assumed.  The 
geometric  configuration,  material  properties  and  computational 
parameters  are  given  in  Fig.  3.4. 

The  evolutions  of  this  stress  wave  propagation  problem  using  IPS2 
element  are  illustrated  in  Fig.  3.4.  One  component  of  the  stress 
tensor  T22  *s  reported  at  various  times.  The  traction  force  h^  »  1000 
is  applied  on  the  boundary  as  a  Heaviside  function  and  this  loading  is 
terminated  at  time  t  -  0.02.  The  Rayleigh  waves  [G raff,  1975]  can  be 
observed  in  Fig.  3.4.  Their  effect  decreases  rapidly  with  depth  and 
their  speed  of  propagation  is  smaller  than  that  of  body  waves. 
Immediately  after  t  -  0.02,  the  finite  element  meshes  are  prescribed  to 
move  with  the  velocity  v^  ■  -0.4/E/p.  At  time  t  ■  0.03,  the 
computations  are  stopped  and  the  stress  distributions  along  lines  y  ■ 
1.5,  2.0,  and  3.0  are  reported  in  Fig.  3.5  for  both  the  elastic  and 
elastic-plastic  cases.  Also  included  in  this  figure  are  the  stresses 
obtained  by  IPS2  and  SRI  elements  without  mesh  motions. 

The  numerical  results  for  the  IPS2  and  SRI  elements  without  mesh 
motion  agree  well  for  both  the  elastic  and  elastic-plastic  cases.  When 
the  mesh  is  moving  with  40%  of  the  elastic  wave  speed,  the  wave  arrival 
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time  agrees  well  while  several  percent  difference  in  wave  amplitudes 
can  oe  observed  as  compared  to  the  fixed  mesh  runs.  This  discrepancy 
is  due  to  the  convective  effects  in  the  momentum  and  constitutive 
equations,  since  these  convective  effects  are  the  only  difference 
between  the  fixed  and  moving  mesh  runs. 
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3.6.3  Plane-Strain  Elastic-Plastic  Beam 


A  two-dimensional  plane-strain  dynamic  elastic-plastic  deformation 
problem  is  used  to  assess  the  proposed  I?S2  element.  The  problem 
statement  and  essential  and  natural  boundary  conditions  are  given  in 
Fig.  3.6.  Small  deformation  is  assumed.  The  entire  beam  is 
discretized  by  25x4  elements.  The  following  material  constants  are 
used:  density  p  -  1;  Young's  modulus  E  -  IQ4;  plastic  tangenc  modulus 
Et  »  0.01E;  Poisson's  ratio  v  *  0.25;  uniaxial  yield  stress  ty  *  300. 
The  Krieg-Key  plasticity  model  [Krieg  and  Key,  1976]  is  employed  and 
isotropic  hardening  Is  assumed. 

The  explicit  predictor-corrector  method  [Hughes  and  Liu,  1978]  is 
employed  with  the  following  computational  parameters:  3  ■  0;  y  •  0.5; 
time  step  size  At  »  0.0075;  and  number  of  time  steps  ■  1000. 

The  time  histories  for  the  tip  displacement  (point  A)  and  stress 
T  at  point  B  are  reported  in  Fig.  3.6.  It  can  be  3een  that  the 
displacement  history  for  the  IPS2  element  is  identical  with  thac  using 
the  SRI  element,  while  several  percent  differences  are  observed  in  the 
stress  history  in  the  plastic  range.  The  system  response  obtained  by 
the  IPS  element  is  also  included  in  Fig.  3.6.  The  maximum  difference 
for  the  displacement  is  approximately  10Z.  To  test  the  sensitivity  of 
the  IPS2  element  to  Irregular  element  shapes,  the  displacement  and 
stress  histories  for  a  fairly  skewed  finite  element  mesh  are  included 
in  Fig.  3.6.  Numerical  solutions  show  that  a  small  amount  of  sciffness 
is  introduced  by  skewing  the  elements. 
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PLANE  STRAIN  ELASTIC -PLASTIC  BEAM 


|  regular  mesh 


B(.5,-1.5)  C(12,0) 

irregular  mesh 


A(25,2) 


Ci(13,— .5) 


p= 1  E=104  i/=0.25 
E/Et=100  Tj=300 
isotropic  hardening 
hy  =  15(l-j*/4) 
jG=.0  7=.5  At=.0075 


(1) IPS2  (reg.)  tg  toe 

(2) IPS2  (irreg.)  g 

(3) SRI  (reg.)  5/5  • 

(4) IPS  (reg.) 


Ic  is  also  worthwhile  to  mention  that  the  computer  time  consumed 
by  the  SRI  element  is  about  3.5  times  of  the  proposed  I?S2  elements 
(using  the  expression  for  internal  force  vector  given  by  Eq.  (3.30)). 

Also  included  in  this  numerical  investigation  are  the  convergence 
studies  for  the  IPS  element.  Fig.  3.7  shows  the  convergence  properties 
for  E/Ej  *  100  and  10,  respectively.  When  the  mesh  is  refined,  the 
responses  converge  to  those  of  the  SRI  element . 


3.6.4  Two  Dimensional  Solid  with  a  Circular  Hole 


A  dynamic  finite  def ormation/ rotation  problem  with  plastic 
hardening  is  considered  here.  As  shown  in  the  proolem  statement  given 
in  Fig.  3.8,  a  plane  strain  square  body  with  a  circular  hole  positioned 
at  the  center  of  the  solid  is  subjected  to  a  uniformly  distributed 
load.  Due  to  symmetry,  only  a  quadrant  of  the  solid  is  modelled  by  360 
elements.  The  dimensions,  material  properties,  essential  and  traction 
boundary  conditions  and  computational  parameters  are  included  in  the 
same  figure.  The  radial  return  procedure  given  in  Krieg  and  Key 
[1976]  is  employed  and  isotropic  hardening  is  assumed.  This  problem  is 
solved  by  degenerating  the  ALE  code  to  the  Lagrangian  description.  In 
other  words,  the  mesh  velocities  are  prescribed  to  be  the  same  as 
material  velocities.  Both  the  SRI  and  IPS2  elements  have  been  tested 
in  this  problem. 

The  dynamic  responses  of  this  system  can  be  observed  from  Fig.  3.9 
in  which  the  mesh  configurations  are  plotted  for  various  times.  The 
x-displacement  for  point  A  and  the  stress  t 22  ^or  Point  3  are  reported 
versus  time  in  Fig.  3.10  for  both  the  SRI  and  IPS2  elements.  The  time 
interval  from  t  -  0  to  t  ~  0.30  can  be  recognized  as  the 
elastic-plastic  loading  period.  The  displacement  and  stress  histories 
for  the  SRI  and  IPS2  elements  agree  very  well  as  can  be  seen  in  Fig. 
3.10.  The  elastic  unloading  occurs  after  t  ~  0.30,  and  a  small  amount 
of  phase  shift  appears  between  these  elements.  However,  both  elements 
exhibit  the  same  shapes  in  the  displacement  and  stress  time  history 
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plots.  After  t  ~  0.75,  the  elastic  loading  occurs  and  the  error 
introduced  in  the  unloading  period  is  retained. 

It  can  be  seen  from  Fig.  3.10  that  the  IPS2  element  is  stiffer  than 
the  SRI  element  (the  same  phenomenon  as  Example  3).  The  computer  time 
consumed  by  the  SRI  element  (on  Harris  800)  is  around  nine  hours  while 
the  IPS2  element  requires  only  2.5  hours  (based  on  the  expression  for 
internal  force  vector  given  by  Eq.  (3.30)).  The  computer  time  ratio 
between  the  SRI  and  the  IPS2  elements  for  this  example  is  approximately 
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3.6.5  Strain  Concentration  with  ALE  Mesh 

A  plane  strain  chick  beam  as  shown  in  Fig.  3.11  is  considered  as  a 
cesc  problem  for  the  ALE  method.  The  dimensions,  material  properties, 
essential  and  natural  boundary  conditions,  and  computational  parameters 
are  given  in  Fig.  3.11.  Constant  shear  loading  is  applied  at  one  end 
of  the  beam  as  a  Heaviside  function.  Small  deformation  is  assumed  and 
only  the  upper  half  of  the  beam  is  modelled  by  the  IPS2  elements 
because  of  the  anti-symmetry  conditions  in  this  example. 

This  example  is  analyzed  by  three  mesh  secups:  (1)  20x4  fixed 
uniform  mesh,  (2)  10x4  fixed  uniform  mesh,  and  (3)  10x4  ALE  mesh.  The 
initial  layout  for  the  10x4  ALE  mexh  is  the  same  as  the  10x4  uniform 
mesh.  As  the  plastic  yielding  effects  are  detected,  the  mesh 
velocities  are  programmed  to  move  according  to  an  ad  hoc  function:  v^  * 
0.15(X-a)^-1.4  |  X-a  |  ,uah*  0  for  0  <  X  <  5,  or  a  *  10  for  5  <  X  ~  10, 
such  that  the  finite  element  mexh  is  concentrated  only  in  the  high 
strain  (or  stress)  regions.  When  the  mesh  size  is  smaller  than  0.5, 
the  mesh  motions  are  terminated  because  of  the  restricton  of  critical 
time  step  for  explicit  time  integration.  The  final  configuration  of 
this  ALE  mesh  is  shown  in  Fig.  3.11.  The  stress  distribution  of 

along  the  line  v  *  1.75  for  these  mesh  setups  are  reported  in  Fig. 
3.12  at  the  instant  of  that  maximum  deflection  occurs  (t  ■  0.9).  As 
can  be. seen  from  these  results,  the  ALE  mesh  provides  fairly  good 
stress  distribution  as  compared  to  that  of  the  20x4  uniform  mesh. 
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